{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Realization of Non-Recursive Filters\n",
    "\n",
    "*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast Convolution\n",
    "\n",
    "The straightforward convolution of two finite-length signals $x[k]$ and $h[k]$ is a numerically complex task. This has led to the development of various techniques with considerably lower complexity. The basic concept of the *fast convolution* is to exploit the correspondence between the convolution and the scalar multiplication in the frequency domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution of Finite-Length Signals\n",
    "\n",
    "The convolution of a causal signal $x_L[k]$ of length $L$ with a causal impulse response $h_N[k]$ of length $N$ is given as\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = x_L[k] * h_N[k] = \\sum_{\\kappa = 0}^{L-1} x_L[\\kappa] \\; h_N[k - \\kappa] = \\sum_{\\kappa = 0}^{N-1} h_N[\\kappa] \\; x_L[k - \\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "where $x_L[k] = 0$ for $k<0 \\wedge k \\geq L$ and $h_N[k] = 0$ for $k<0 \\wedge k \\geq N$. The resulting signal $y[k]$ is of finite length $M = N+L-1$. The computation of $y[k]$ for $k=0,1, \\dots, M-1$ requires $M \\cdot N$ multiplications and $M \\cdot (N-1)$ additions. The computational complexity of the convolution is consequently [in the order of](https://en.wikipedia.org/wiki/Big_O_notation) $\\mathcal{O}(M \\cdot N)$. Discrete-time Fourier transformation (DTFT) of above relation yields\n",
    "\n",
    "\\begin{equation}\n",
    "Y(e^{j \\Omega}) = X_L(e^{j \\Omega}) \\cdot H_N(e^{j \\Omega})\n",
    "\\end{equation}\n",
    "\n",
    "Discarding the effort of transformation, the computationally complex convolution is replaced by a scalar multiplication with respect to the frequency $\\Omega$. However, $\\Omega$ is a continuous frequency variable which limits the numerical evaluation of this scalar multiplication. In practice, the DTFT is replaced by the discrete Fourier transformation (DFT). Two aspects have to be considered before a straightforward application of the DFT\n",
    "\n",
    "1. The DFTs $X_L[\\mu]$ and $H_N[\\mu]$ are of length $L$ and $N$ respectively and cannot be multiplied straightforwardly\n",
    "    \n",
    "2. For $N = L$, the multiplication of the two spectra $X_L[\\mu]$ and $H_L[\\mu]$ would result in the [periodic/circular convolution](https://en.wikipedia.org/wiki/Circular_convolution) $x_L[k] \\circledast h_L[k]$ due to the periodicity of the DFT. Since we aim at realizing the linear convolution $x_L[k] * h_N[k]$ with the DFT, special care has to be taken to avoid cyclic effects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Convolution by Periodic Convolution\n",
    "\n",
    "The periodic convolution of the two signals $x_L[k]$ and $h_N[k]$ is defined as\n",
    "\n",
    "\\begin{equation}\n",
    "x_L[k] \\circledast h_N[k] = \\sum_{\\kappa=0}^{M-1} \\tilde{x}_M[k - \\kappa] \\; \\tilde{h}_M[\\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "where the periodic continuations $\\tilde{x}_M[k]$ of $x_L[k]$ and $\\tilde{h}_M[k]$ of $h_N[k]$ with period $M$ are given as\n",
    "\n",
    "\\begin{align}\n",
    "\\tilde{x}_M[k] &= \\sum_{m = -\\infty}^{\\infty} x_L[m \\cdot M + k] \\\\\n",
    "\\tilde{h}_M[k] &= \\sum_{m = -\\infty}^{\\infty} h_N[m \\cdot M + k]\n",
    "\\end{align}\n",
    "\n",
    "The result of the circular convolution has a periodicity of $M$.\n",
    "\n",
    "To compute the linear convolution by the periodic convolution one has to take care that the result of the linear convolution fits into one period of the periodic convolution. Hence, the periodicity has to be chosen as $M \\geq N+L-1$. This can be achieved by zero-padding of $x_L[k]$ and $h_N[k]$ to a total length of $M$\n",
    "\n",
    "\\begin{align}\n",
    "x_M[k] &= \\begin{cases}\n",
    "x_L[k] & \\mathrm{for} \\; k=0, 1, \\dots, L-1 \\\\\n",
    "0 & \\mathrm{for} \\; k=L, L+1, \\dots, M-1\n",
    "\\end{cases}\n",
    "\\\\\n",
    "h_M[k] &= \\begin{cases}\n",
    "h_N[k] & \\mathrm{for} \\; k=0, 1, \\dots, N-1 \\\\\n",
    "0 & \\mathrm{for} \\; k=N, N+1, \\dots, M-1\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "This results in the desired equality of linear and periodic convolution\n",
    "\n",
    "\\begin{equation}\n",
    "x_L[k] * h_N[k] = x_M[k] \\circledast h_M[k]\n",
    "\\end{equation}\n",
    "\n",
    "for $k = 0,1,\\dots, M-1$ with $M = N+L-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Linear by periodic convolution\n",
    "\n",
    "The following example computes the linear, periodic and linear by periodic convolution of a rectangular signal $x[k] = \\text{rect}_L[k]$ of length $L$ with a triangular signal $h[k] = \\Lambda_N[k]$ of length $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAADkCAYAAAAsAtjDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFH5JREFUeJzt3X/0nnV93/Hnq0lQNoqUxnowiSTdUkvGD+nJkB7dKVV2BNoD1tEOVlp/rehWXHdKtbBuWFndVjkrHlemZq2COPmhU5pxsmYKeFatIsFoIGSpkdUSoOIPQFv5keB7f9xX4ObLnXzvL36/1/dz534+zvme3Nd1ve/P9b64Dt+8cv1MVSFJkqR2/NBiNyBJkqSnM6BJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJImUpLtSU7pYT1/meTUZ/G9SvK3Sd45Y/49SU4cUX9zkkeTfOYH6VfSwcGAJqlZSV6e5M+TPJzk20k+m+QfAlTVP6iqTy9yi7M5oap+Z99Ekh8BjgJ2zCysqlcAb+6xN0kNW7rYDUjSKEkOB24E/gVwPXAI8I+Axxazrx/QccCuqnp0sRuR1DaPoElq1U8AVNU1VfVEVT1SVf+7qrbB0089JvmpJFuTfDfJR5Ncl+T39g3U1f5Wkm3d0bjrkjx3aPlFSb7aff+uJL8wToNJ3pXkE0PTlyW5Kcmy/XzleODOrvbvJPlIko8nOWyu/3EkHdwMaJJa9RfAE0muSnJ6d3rwGZIcAnwCuBI4ErgGGBWwfgk4DVjDICi9bmjZVxkcnXse8A7gw0mOGqPH3wd+NslLkry5G/81VbVnP/XHA3ckWQN8BtgJ/JOq+psx1iVpihjQJDWpqr4DvBwo4L8B30iyMckLZpSezOByjfdU1Z6q+jjwhRFDvqeq7quqbwP/E3jJ0Lo+2i37flVdB3wFOGmMHr8FvBv4EHAxcEZVPXyArxzH4Bq0m4F3VNU7qqpmW4+k6WNAk9SsqtpRVa+rqpXAscALGQSiYS8E7p0RdO4ZMdxfD33+HvDkacUkv5rkS0keSvJQt67lY7a5lUHwuriqRq133zrSjfsLwPuq6k/GHF/SFDKgSZoIVfV/GZzGPHbGovuBFV0A2mfVuOMmOZrBEboLgB+tqiMYXCeWA35x8N3jgPcCVwFvmKV8TffnqcCFSdaP26Ok6WNAk9SkJD+Z5MIkK7vpVcC5wOdnlH4OeAK4IMnSJGcxxunJIX+XwWnUb3TreT3PDIGj+lvB4FTpm4F/CRw3y3PZjge2VdUdwPnAJ8a8zk3SFDKgSWrVd4GXArcm+VsGwexO4MLhoqp6HHgN8EbgIeA8Bo/nGOtxHFV1F/CfGQS9rzM4XfnZA32newTIJuAPqmpjVX0PuAx45wG+dhywrVvnDcAG4Ibhu0klaZ94faqkg02SWxlc5/XBRezhUQYh8T1V9e/GqP8kgxsevlBVr1zo/iS1zYAmaeIl+RkGj6z4JvDLwPuAH6+q+xe1MUl6lnyTgKSDwYsZvG3gMAbPNDvbcCZpknkETZIkqTHeJCBJktQYA5okSVJjJv4atOXLl9fq1asXuw1JkqRZ3X777d+squfPVjfxAW316tVs2bJlsduQJEmaVZKvjVPnKU5JkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SZKkxhjQJEmSGmNAkyRJaowBTZIkqTEGNEmSpMb0FtCSfCDJA0nu3M/yJHlPkl1JtiX5qb56kyRJakmfr3q6EvhD4EP7WX46sLb7eSnw3u7PA7rj3od52X+6mbe+6sW8+sQV+627Yeu9XLZ5J/c99AgvPOLQA9bPpXah6yd1bHuZ7u1sqZdp2c6Wepnr2JKeKVXV38qS1cCNVXXsiGXvBz5dVdd00zuBU6rq/gON+Zyj1tZRr303hy5bwn98zXEjfwncsPVeLv74HTyy54kn5+2vfi61C10/qWPby3RvZ0u9TMt2ttTLXMeWpk2S26tq/Wx1LV2DtgK4Z2h6dzdvLI/seYLLNu8cueyyzTuf9sviQPVzqV3o+kkd216meztb6mVatrOlXuY6tqTRWgpoGTFv5OG9JOcn2ZJky/D8+x56ZOTAc5k/H2NMQi/Tsp0t9TIt29lSL9OynS31MtcxJI3WUkDbDawaml4J3DeqsKo2VNX6mYcIX3jEoSMHnsv8+RhjEnqZlu1sqZdp2c6WepmW7Wypl7mOIWm0lgLaRuBXu7s5TwYenu36s2GHLlvCW1/14pHL3vqqF3PosiVj1c+ldqHrJ3Vse5nu7Wypl2nZzpZ6mevYkkbr7S7OJNcApwDLk+wG3g4sA6iq9wGbgDOAXcD3gNePO/aKWe4S2jf/bR/bxuNPfP+A9XOpXej6SR3bXqZ7O1vqZVq2s6Ve5jq2pNF6vYtzIRx59DH17a/tGKv2n77/cwBc96afntfaha6f1LHtpf+x7aX/se1lfsaWpsUk3sUpSZIkDGiSJEnNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjDGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmN6DWhJTkuyM8muJBeNWP6iJLck2ZpkW5Iz+uxPkiSpBb0FtCRLgCuA04F1wLlJ1s0o+7fA9VV1InAO8F/76k+SJKkVfR5BOwnYVVV3V9XjwLXAWTNqCji8+/w84L4e+5MkSWpCnwFtBXDP0PTubt6w3wXOS7Ib2AS8ZdRASc5PsiXJlj179ixEr5IkSYumz4CWEfNqxvS5wJVVtRI4A7g6yTN6rKoNVbW+qtYvW7ZsAVqVJElaPH0GtN3AqqHplTzzFOYbgesBqupzwHOB5b10J0mS1Ig+A9ptwNoka5IcwuAmgI0zav4KeCVAkmMYBLRv9NijJEnSoustoFXVXuACYDOwg8HdmtuTXJrkzK7sQuDXknwZuAZ4XVXNPA0qSZJ0UFva58qqahODi/+H510y9Pku4GV99iRJktQa3yQgSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjDGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDWm14CW5LQkO5PsSnLRfmp+KcldSbYn+Uif/UmSJLVgaV8rSrIEuAL4x8Bu4LYkG6vqrqGatcDFwMuq6sEkP9ZXf5IkSa3o8wjaScCuqrq7qh4HrgXOmlHza8AVVfUgQFU90GN/kiRJTZj1CFqSI8cY5/tV9dAsNSuAe4amdwMvnVHzE906PwssAX63qv50jPVLkiQdNMY5xXlf95MD1CwBXjTLOKO+XyP6WQucAqwE/izJsTPDX5LzgfMBDjvq782yWkmSpMkyTkDbUVUnHqggydYxxtkNrBqaXskg+M2s+XxV7QH+X5KdDALbbcNFVbUB2ABw5NHHzAx5kiRJE22ca9B+GiDJ781c0F34/2TNLG4D1iZZk+QQ4Bxg44yaG4Cf7cZezuCU591jjC1JknTQmDWgVdWj3ccVSf7ZvvndHZafmlFzoHH2AhcAm4EdwPVVtT3JpUnO7Mo2A99KchdwC/DWqvrWXDZIkiRp0s3lMRtvAjYn2cXg2rEPAr89l5VV1SZg04x5lwx9LuA3ux9JkqSpNM5dnB8CvghsBX4d+AiwF3h1Ve1a2PYkSZKmzzjXoF3V1b2BQThbDTwInJfk7IVrTZIkaTrNegStqm4Cbto3nWQpsA44ATgZ+NiCdSdJkjSF5vyqp+5i/23dz9Xz3pEkSdKUm/UUZ5IvzkeNJEmSxjPOEbRjkmw7wPIAz5unfiRJkqbeOAHtJ8eoeeIHbUSSJEkD49wk8DWAJJ8CLqyqLy94V5IkSVNsnMds7PM24PIkH0xy1EI1JEmSNO3GDmhV9cWqegVwI/CnSd6e5NCFa02SJGk6zeUIGkkC7ATeC7wF+EqSX1mIxiRJkqbV2AEtyWeAe4HLgRXA64BTgJOSbFiI5iRJkqbRXB5U+2Zge/dC82FvSbJjHnuSJEmaamMHtKq68wCLf24eepEkSRJzvAZtf6rq7vkYR5IkSfMU0CRJkjR/DGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjek1oCU5LcnOJLuSXHSAurOTVJL1ffYnSZLUgt4CWpIlwBXA6cA64Nwk60bU/TDwr4Bb++pNkiSpJX0eQTsJ2FVVd1fV48C1wFkj6v498C7g0R57kyRJakafAW0FcM/Q9O5u3pOSnAisqqobDzRQkvOTbEmyZc+ePfPfqSRJ0iLqM6BlxLwnX7ye5IeAy4ELZxuoqjZU1fqqWr9s2bJ5bFGSJGnx9RnQdgOrhqZXAvcNTf8wcCzw6SR/CZwMbPRGAUmSNG36DGi3AWuTrElyCHAOsHHfwqp6uKqWV9XqqloNfB44s6q29NijJEnSoustoFXVXuACYDOwA7i+qrYnuTTJmX31IUmS1Lqlfa6sqjYBm2bMu2Q/taf00ZMkSVJrfJOAJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjDGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjeg1oSU5LsjPJriQXjVj+m0nuSrItyU1Jju6zP0mSpBb0FtCSLAGuAE4H1gHnJlk3o2wrsL6qjgc+Bryrr/4kSZJa0ecRtJOAXVV1d1U9DlwLnDVcUFW3VNX3usnPAyt77E+SJKkJfQa0FcA9Q9O7u3n780bgf41akOT8JFuSbNmzZ888tihJkrT4+gxoGTGvRhYm5wHrgctGLa+qDVW1vqrWL1u2bB5blCRJWnxLe1zXbmDV0PRK4L6ZRUlOBX4H+Jmqeqyn3iRJkprR5xG024C1SdYkOQQ4B9g4XJDkROD9wJlV9UCPvUmSJDWjt4BWVXuBC4DNwA7g+qranuTSJGd2ZZcBhwEfTfKlJBv3M5wkSdJBq89TnFTVJmDTjHmXDH0+tc9+JEmSWuSbBCRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SZKkxhjQJEmSGmNAkyRJaowBTZIkqTEGNEmSpMYY0CRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SZKkxvQa0JKclmRnkl1JLhqx/DlJruuW35pkdZ/9SZIktaC3gJZkCXAFcDqwDjg3yboZZW8EHqyqvw9cDvx+X/1JkiS1os8jaCcBu6rq7qp6HLgWOGtGzVnAVd3njwGvTJIee5QkSVp0qap+VpScDZxWVf+8m/4V4KVVdcFQzZ1dze5u+qtdzTf3N+7aw4+oT5115lg93HX/dwBYd9Th81q70PWTOra99D+2vfQ/tr3sv/avn7+K11/zX8YaW5oWSW6vqvWz1S3to5nOqCNhM9PhODUkOR84v5t8bPWHr77zB+xN7VgO7DeQayK5Tw8uc9qfb7j2DxewFc0D///s39HjFPUZ0HYDq4amVwL37admd5KlwPOAb88cqKo2ABsAkmwZJ4lqMrg/Dz7u04OL+/Pg4v5sV5/XoN0GrE2yJskhwDnAxhk1G4HXdp/PBm6uvs7BSpIkNaK3I2hVtTfJBcBmYAnwgaranuRSYEtVbQT+GLg6yS4GR87O6as/SZKkVvR5ipOq2gRsmjHvkqHPjwK/OMdhN8xDa2qH+/Pg4z49uLg/Dy7uz0b1dhenJEmSxuOrniRJkhoz0QFttldHqW1JPpDkge75d/vmHZnkk0m+0v35I4vZo8aXZFWSW5LsSLI9yW90892nEyjJc5N8IcmXu/35jm7+mu5VfF/pXs13yGL3qvElWZJka5Ibu2n3Z6MmNqCN+eoote1K4LQZ8y4CbqqqtcBN3bQmw17gwqo6BjgZ+PXu/0n36WR6DHhFVZ0AvAQ4LcnJDF7Bd3m3Px9k8Io+TY7fAHYMTbs/GzWxAY3xXh2lhlXV/+GZz7kbft3XVcCre21Kz1pV3V9VX+w+f5fBXwIrcJ9OpBr4m25yWfdTwCsYvIoP3J8TJclK4OeAP+qmg/uzWZMc0FYA9wxN7+7mabK9oKruh8Ff+MCPLXI/ehaSrAZOBG7FfTqxutNhXwIeAD4JfBV4qKr2diX+3p0s7wbeBny/m/5R3J/NmuSANtZroST1K8lhwP8A/nVVfWex+9GzV1VPVNVLGLz55STgmFFl/XalZyPJzwMPVNXtw7NHlLo/G9Hrc9Dm2TivjtLk+XqSo6rq/iRHMfiXuyZEkmUMwtl/r6qPd7PdpxOuqh5K8mkG1xYekWRpd9TF37uT42XAmUnOAJ4LHM7giJr7s1GTfARtnFdHafIMv+7rtcCfLGIvmoPuepY/BnZU1R8MLXKfTqAkz09yRPf5UOBUBtcV3sLgVXzg/pwYVXVxVa2sqtUM/r68uap+Gfdnsyb6QbXdvwTezVOvjnrnIrekOUhyDXAKsBz4OvB24AbgeuBFwF8Bv1hVM28kUIOSvBz4M+AOnrrG5d8wuA7NfTphkhzP4KLxJQz+MX99VV2a5McZ3JR1JLAVOK+qHlu8TjVXSU4Bfquqft792a6JDmiSJEkHo0k+xSlJknRQMqBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiQNSXJqkqsXuw9J082AJklPdwKDJ6pL0qIxoEnS050AbE3ynCRXJvkP3XtGJak3Sxe7AUlqzAnAA8Bm4I+q6sOL3I+kKeS7OCWpk2QZ8E3ga8Cbqupzi9ySpCnlKU5Jeso64DZgL/DEIvciaYoZ0CTpKScAfw6cA3wwyQsWuR9JU8qAJklPOQG4s6r+Avht4PrutKck9cpr0CRJkhrjETRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SZKkxhjQJEmSGmNAkyRJaowBTZIkqTH/HyseHO3/mGY6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAADiCAYAAABwdKKfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF09JREFUeJzt3X+0XWV54PHv4yWUixRSITrkggZaJ2NW+ZGuLAWpHRdlVqgwNWUxI1q11nGga4ljpxhKOp1V2w5Km45a66wuI4pUbBomzWQcYRmtwBqtDhi4mCCZq4iA3FAJ6lVa75Dk8swfZwdv4s3N2TfnvOfsc76ftc7KOfu8e+/nZGe/75O93/2+kZlIkiSpjOf1OgBJkqRhYvIlSZJUkMmXJElSQSZfkiRJBZl8SZIkFWTyJUmSVJDJl6RaIuJVETHR6zj6UUTcFRFvO4r1/zEizuxkTJL6j8mXpDlFxCMRcdGhyzPzC5m5vBcxDZK5ErXMPCEzH+5VTJLKMPmS1AgRcUyvY5CkTjD5klRLRLw6Ih6f9fmRiHhXROyIiB9ExKaIOG7W95dGxP0RMRURX4qIs2d9d11EfDMino6IByPi12Z995aI+PuIeH9EfA949xyxjETE783axr0RcXr13Ssj4itVTF+JiFfOWu+uiPjjavtPR8RnI+KU6rvPRMTVh+znqxFx2ZG2e8g6746IW2Z9XhYRGRHHRMT1wKuAD1W3Gj9UlcmI+Lnq/UkR8VcRsSciHo2I34+I5836u/liRPxZRHw/Ir4VEb/SzvGT1HsmX5I64d8CFwNnAGcDbwGIiF8APgZcBZwMfBj4VET8VLXeN2klIScBfwjcEhGnztruK4CHgRcC18+x398BXg+8BjgReCvwo4h4AXAb8MFqv+8DbouIk2et+wbgN6ttHwu8q1r+19U2qX7DCuAl1frtbPeIMvM/AV8Arq5uNV49R7G/oPX3cibwL4E3V/Ee8ApgAjgF+FPgoxERdeKQ1BsmX5I64YOZuTszvwf8L+Dcavm/Bz6cmXdn5kxm3gw8A5wHkJn/vVrv2czcBHwDePms7e7OzL/IzP2ZOT3Hft8G/H5mTmTLVzPzu8AlwDcy8xPVuhuB/wv861nr3pSZX6+2e+usmP8HcG5EvKT6/OvAlsx8ps3tHrWIGAFeB6zLzKcz8xHgvwJvmlXs0cz8SGbOADcDpwIv6mQckrrD5EtSJ/zDrPc/Ak6o3r8EuKa65TgVEVPA6cBSgIh486xbklPAz9O6knPAt4+w39NpXT071FLg0UOWPQqMHSnmzHya1tWtK6rvrgA+WWO7nXAKratxs/d12Pgz80fV2xOQ1PdMviR107eB6zNz8azX8Zm5sbqy9BHgauDkzFwMPADMvnWWbWz/Z+dYvptW4jfbi4HJNuPeCLw+Is4HRoE7F7DdfwKOn/X5nx3y/Xy/7Slg3yH7qhO/pD5m8iVpPosi4rhZr7pPHH4E+K2IeEW0PD8iLomInwaeTysB2QMQEb9J68pXHTcCfxwRL622f3bV/+p24J9HxBuqDu6vA1YAn25zu7fTSnz+CNiUmc/OWt7udu8HfikiXhwRJwHrDvn+O7T6c/2E6lbircD1EfHTVaL6O8Atc5WX1CwmX5LmczswPev17jorZ+Z2Wv2+PgR8H3iIqjN+Zj5Iqx/Tl2klImcBf18zvvfRSlI+C/wQ+CgwWvX7uhS4BvgucC1waWY+1WbczwBbgItodcA/sLzt7Wbm54BNwA7gXn4yQftz4PLqacUPzhHGO2hdPXsY+GIVx8faiV9Sf4vMI13VlyRJUqd45UuSJKkgky9JkqSCTL4kSZIKMvmSJEkqyORLkiSpoLpj9hR1yimn5LJly3odhiRJ0hHde++9T2XmkiOV6+vka9myZWzfvr3XYUiSJB1RRBw6/dicvO0oSZJUkMmXJElSQSZfkiRJBZl8SZIkFWTyJUmSVJDJlyRJUkEmX5IkSQWZfEmSJBVk8iVJklRQ0eQrIv5jRHwtIh6IiI0RcVzJ/UuSJPVasemFImIM+A/AisycjohbgSuAj5eKQeqmreOTrN82we6paZYuHmXt6uWsWTnWkfJN3baxDPfvlDS30nM7HgOMRsQ+4Hhgd+H9S22r24Ct27KT6X0zAExOTbNuy06AOdepU76p2zaW4f6ds9cxWZMOVuy2Y2ZOAn8GPAY8AfwgMz9bav9SHQcamcmpaZIfNzJbxyfnLL9+28RzDdIB0/tmWL9t4qjLN3XbxjLcvxPqn0fSsCiWfEXEzwCvBc4AlgLPj4g3zlHuyojYHhHb9+zZUyo86SB1G5ndU9NdW97UbRvLcP9OqH8eScOiZIf7i4BvZeaezNwHbAFeeWihzNyQmasyc9WSJUsKhqdhsHV8kgtuuIMzrruNC26447D/A6/byCxdPNq15U3dtrEM9++E+udRu+en1HQlk6/HgPMi4viICOCXgV0F968hV+cWSN1GZu3q5YwuGjlo2eiiEdauXn7U5Zu6bWMZ7t8J9c4jb1FqmJTs83U3sBm4D9hZ7XtDqf1LdW6B1G1k1qwc472XncWxI61TamzxKO+97KzDdiyuU76p2zaW4f6dUO888halhklkZq9jOKxVq1bl9u3bex2GBsQZ193GXP/aA/jWDZf8xPKt45Ncu3kHe2eeZazNp7Re9+EvA7DpqvPbiqlO+aZu21jKb7ufYmn3PKp7fkr9KCLuzcxVRypXeqgJqWeWLh5lco6+Joe7NbJm5Rgb73kMaL9RknSwds+juuen1GROL6ShUfdWoqRyPD81TLzypUarM4DjgeV1byVK6r6FnJ8O4KqmMvlSYy1ktG1vJUr9q875uZDzX+oX3nZUY/l0lDS8PP/VZCZfaqy6AzhKGhye/2oyky81Vt2BUCUNDs9/NZnJlxrLp6Ok4eX5ryazw70ay6cXpeHl+a8mM/lSo/n0ojS8PP/VVCZf6iuO2yOpG6xb1E9MvtQ3HLdHUjdYt6jf2OFefcNxeyR1g3WL+o3Jl/qG4/ZI6gbrFvUbky/1DcftkdQN1i3qNyZf6huO2yOpG6xb1G/scK++4bg9krrBukX9xuRLfcVxeyR1g3WL+om3HSVJkgoy+ZIkSSrI5EuSJKkg+3yp65zWQ1LTWG+pm0y+1FVO6yGpaay31G3edlRXOa2HpKax3lK3mXypq5zWQ1LTWG+p20y+1FVO6yGpaay31G0mX+oqp/WQ1DTWW+o2O9yrq5zWQ1LTWG+p24omXxGxGLgR+Hkggbdm5pdLxqDynNZDUtNYb6mbSl/5+nPgM5l5eUQcCxxfeP+SJEk9VSz5iogTgV8C3gKQmXuBvaX2L0mS1A9Kdrg/E9gD3BQR4xFxY0Q8v+D+JUmSeq5k8nUM8AvAX2bmSuCfgOsOLRQRV0bE9ojYvmfPnoLhSZIkdV/JPl+PA49n5t3V583MkXxl5gZgA8CqVauyXHhql3OeSdLBrBdVR7HkKzP/ISK+HRHLM3MC+GXgwVL7V2c455kkHcx6UXWVHmT1HcAnI2IHcC7wnsL711FyzjNJOpj1ouoqOtREZt4PrCq5T3WWc55J0sGsF1WX0wupFuc8k6SDWS+qLpMv1eKcZ5J0MOtF1eXcjqrFOc8k6WDWi6rL5Eu1OeeZJB3MelF1eNtRkiSpIJMvSZKkgky+JEmSCjL5kiRJKsjkS5IkqSCTL0mSpIJMviRJkgpynC8BsHV8kvXbJtg9Nc1SBwiUpK6xvpXJl9g6Psm6LTuZ3jcDwOTUNOu27ASwQpCkDrK+FXjbUcD6bRPPVQQHTO+bYf22iR5FJEmDyfpWYPIlYPfUdK3lkqSFsb4VmHwJWLp4tNZySdLCWN8K2ki+IuIFbbwWlwhW3bF29XJGF40ctGx00QhrVy/vUUSSNJisbwXtdbjfXb1injIjwIs7EpGKO9DJ89rNO9g78yxjPn0jSV1hfStoL/nalZkr5ysQEeMdikc9smblGBvveQyATVed3+NoJGlwWd+qnT5f5wNExH859IuIGJldRpIkSfM7YvKVmf+vejsWEW84sDwiXgj83SFlJEmSNI86g6xeBWyLiIeABG4CfrcrUUmSJA2oIyZfEfFXwH3AOPB24K+B/cCazHyou+FJkiQNlnb6fN1clXsrrcRrGfB94I0RcXn3QpMkSRo8R7zylZmfBz5/4HNEHAOsAM4BzgM2dy06SZKkAVN7Yu3M3A/sqF6f6HhEkiRJA6ydEe7v60QZSZIktXfl62URsWOe7wM4qUPxSJIkDbR2kq9/0UaZmXZ3WA3Muh2YzMxL211PkiRpELTT4f5RgIj4O+CazPzqUe7zncAu4MSj3I7msXV8kvXbJtg9Nc1S5w6TpMayPh887Qw1ccC1wPsj4qaIOHUhO4uI04BLgBsXsr7as3V8knVbdjI5NU0Ck1PTrNuyk63jk70OTZJUg/X5YGo7+crM+zLzQuDTwGci4g8iYrTm/j5AK4l7tuZ6qmH9tgmm9x18J3h63wzrt030KCJJ0kJYnw+mOle+iIgAJoC/BN4BfCMi3tTmupcCT2bmvUcod2VEbI+I7Xv27KkTniq7p6ZrLZck9Sfr88HUdvIVEV8EJoH3A2PAW4BXAy+PiA1tbOIC4Fcj4hHgb4ALI+KWQwtl5obMXJWZq5YsWdJueJpl6eK5L0gebrkkqT9Znw+mOle+fgsYy8x/lZn/OTM/nZkPZeY7gFcdaeXMXJeZp2XmMuAK4I7MfOPCwtZ81q5ezuiikYOWjS4aYe3q5T2KSJK0ENbng6ntEe4z84F5vr6kA7GoQw48BXPt5h3snXmWMZ+OkaRGsj4fTLWnF5pLZj5cs/xdwF2d2LfmtmblGBvveQyATVed3+NoJEkLZX0+eGp1uJckSdLRMfmSJEkqyORLkiSpIJMvSZKkgky+JEmSCjL5kiRJKsjkS5IkqSCTL0mSpIJMviRJkgrqyAj36r6t45Os3zbB7qlpljq9hCTpMGwv+p/JVwNsHZ9k3ZadTO+bAWByapp1W3YCeEJJkp5je9EM3nZsgPXbJp47kQ6Y3jfD+m0TPYpIktSPbC+aweSrAXZPTddaLkkaTrYXzWDy1QBLF4/WWi5JGk62F81g8tUAa1cvZ3TRyEHLRheNsHb18h5FJEnqR7YXzWCH+wY40Eny2s072DvzLGM+vSJJmoPtRTOYfDXEmpVjbLznMQA2XXV+j6ORJPUr24v+521HSZKkgky+JEmSCjL5kiRJKsjkS5IkqSCTL0mSpIJMviRJkgoy+ZIkSSrIcb56aOv4JOu3TbB7apqlDoQnSSrMdqg3TL56ZOv4JOu27Hxu9vnJqWnWbdkJ4D98SVLX2Q71jrcde2T9tonn/sEfML1vhvXbJnoUkSRpmNgO9Y7JV4/snpqutVySpE6yHeqdYslXRJweEXdGxK6I+FpEvLPUvvvR0sWjtZZLktRJtkO9U/LK137gmsx8GXAe8PaIWFFw/31l7erljC4aOWjZ6KIR1q5e3qOIJEnDxHaod4p1uM/MJ4AnqvdPR8QuYAx4sFQM/eRAZ8ZrN+9g78yzjPmUiSSpINuh3unJ044RsQxYCdzdi/33izUrx9h4z2MAbLrq/B5HI0kaNrZDvVG8w31EnAD8LfDbmfnDOb6/MiK2R8T2PXv2lA5PkiSpq4omXxGxiFbi9cnM3DJXmczckJmrMnPVkiVLSoYnSZLUdSWfdgzgo8CuzHxfqf1KkiT1k5JXvi4A3gRcGBH3V6/XFNy/JElSz5V82vGLQJTanyRJUj9ybscOcoJSSdIgs53rDJOvDnGCUknSILOd6xznduwQJyiVJA0y27nOMfnqECcolSQNMtu5zjH56hAnKJUkDTLbuc4x+eoQJyiVJA0y27nOscN9hzhBqSRpkNnOdY7JVwc5QakkaZDZznWGtx0lSZIKMvmSJEkqyNuO83AkX0mSFsY29PBMvg7DkXwlSVoY29D5edvxMBzJV5KkhbENnZ/J12E4kq8kSQtjGzo/k6/DcCRfSZIWxjZ0fiZfh+FIvpIkLYxt6PzscH8YjuQrSdLC2IbOz+RrHo7kK0nSwtiGHt7QJV+OOyJJUv8ZpvZ5qJIvxx2RJKn/DFv7PFQd7h13RJKk/jNs7fNQJV+OOyJJUv8ZtvZ5qJIvxx2RJKn/DFv7PFTJl+OOSJLUf4atfW58h/s6T0c47ogkSf1nIe1zk5+ObHTytZCnIxx3RJKk/lOnfW7605GNvu04bE9HSJKk5rf/jU6+hu3pCEmS1Pz2v2jyFREXR8RERDwUEdcdqfzOyR9wwQ13sHV8cs7vh+3pCEmSVL/93zo+yQU33MEZ1902b15RSrHkKyJGgP8G/AqwAnh9RKw40noH7uPO9Rc1bE9HSJKkeu3/gf5hk1PTJPPnFaWUvPL1cuChzHw4M/cCfwO8tp0VD3cfd83KMd572VkcO9L6GWOLR3nvZWc1orOdJElamDrtfz/2D4vMLLOjiMuBizPzbdXnNwGvyMyrD7fOz55wUr7nnF987vN5Z548Z7kHn/ghACtOPbGtWOqU7+a2jaX8to2l/LaNpfy2jaX8to2l/LbbLf9/Hv7uc+8fPmmMD5/duu4TwLduuKSt/bQrIu7NzFVHKldyqImYY9lPZH4RcSVwJcDzRk/kzQ9/vVVwZv/efV96ZGdXI1S3nQI81esg1DEez8Hi8Rw8HlNg0ZJlZ8XIMce2Pn0dxu8EWnlF/Mmlnc4rXtJOoZLJ1+PA6bM+nwbsPrRQZm4ANgBExPZnfvSDI2aQaoaI2N7O/wjUDB7PweLxHDwe0/5Vss/XV4CXRsQZEXEscAXwqYL7lyRJ6rliV74yc39EXA1sA0aAj2Xm10rtX5IkqR8UnV4oM28Hbq+xyoZuxaKe8HgOFo/nYPF4Dh6PaZ8q9rSjJEmSGj69kCRJUtP0ZfJVdxoi9Z+I+FhEPBkRD8xa9oKI+FxEfKP682d6GaPaFxGnR8SdEbErIr4WEe+slntMGygijouIeyLiq9Xx/MNq+RkRcXd1PDdVD0epISJiJCLGI+LT1WePZ5/qu+RrodMQqe98HLj4kGXXAZ/PzJcCn68+qxn2A9dk5suA84C3V+elx7SZngEuzMxzgHOBiyPiPOBPgPdXx/P7wL/rYYyq753ArlmfPZ59qu+SL45iGiL1j8z838D3Dln8WuDm6v3NwJqiQWnBMvOJzLyvev80rQp+DI9pI2XLP1YfF1WvBC4ENlfLPZ4NEhGnAZcAN1afA49n3+rH5GsM+Pasz49Xy9R8L8rMJ6DVmAMv7HE8WoCIWAasBO7GY9pY1S2q+4Engc8B3wSmMnN/VcS6t1k+AFwLPFt9PhmPZ9/qx+SrrWmIJJUXEScAfwv8dmb+sNfxaOEycyYzz6U128jLgZfNVaxsVFqIiLgUeDIz7529eI6iHs8+UXScrza1NQ2RGuk7EXFqZj4REafS+h+3GiIiFtFKvD6ZmVuqxR7ThsvMqYi4i1ZfvsURcUx1tcS6tzkuAH41Il4DHAecSOtKmMezT/XjlS+nIRpcnwJ+o3r/G8D/7GEsqqHqP/JRYFdmvm/WVx7TBoqIJRGxuHo/ClxEqx/fncDlVTGPZ0Nk5rrMPC0zl9FqM+/IzF/H49m3+nKQ1Sp7/wA/nobo+h6HpJoiYiPwauAU4DvAHwBbgVuBFwOPAf8mMw/tlK8+FBG/CHwB2MmP+5T8Hq1+Xx7ThomIs2l1wB6h9Z/wWzPzjyLiTFoPOb0AGAfemJnP9C5S1RURrwbelZmXejz7V18mX5IkSYOqH287SpIkDSyTL0mSpIJMviRJkgoy+ZIkSSrI5EuSJKkgky9JkqSCTL4kSZIKMvmSNDQi4qKI+ESv45A03Ey+JA2Tc2iN9C1JPWPyJWmYnAOMR8RPRcTHI+I91byVklTMMb0OQJIKOgd4EtgG3JiZt/Q4HklDyLkdJQ2FiFgEPAU8ClyVmV/ucUiShpS3HSUNixXAV4D9wEyPY5E0xEy+JA2Lc4AvAVcAN0XEi3ocj6QhZfIlaVicAzyQmV8Hfhe4tboVKUlF2edLkiSpIK98SZIkFWTyJUmSVJDJlyRJUkEmX5IkSQWZfEmSJBVk8iVJklSQyZckSVJBJl+SJEkF/X+enPC1Mbba3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAADiCAYAAAAYhyCnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG7lJREFUeJzt3X+8ZXVd7/HXu5mDkoCAIPIbU0Kg26DNRZS0yQyREOxeLNAUfzVS2iPvtdSsC2jWozKlR2rRJD8UFX+kIikIqAlyM3X4DQKCXJBhRkZ+CSbqDH7uH2tNbvbsM2fPcM7ee53zej4e+3HW+q7vWvuz9/fAec/6mapCkiRJ3fAz4y5AkiRJwzO8SZIkdYjhTZIkqUMMb5IkSR1ieJMkSeoQw5skSVKHGN6kDkjykiQXbuG6Jyf5YDu9V5LvJ1k0uxWOVpIzk7z9Eax/XZJls1jSrL5/ki8lefUIS+p972cluXEL1315kktnuyZJD2d4k+ZIkluTPNiGpTuTnJFkmy3ZVlV9qKoOe6Q1VdW3q2qbqnrokW6rKwYFvao6sKq+NKaSHvb+veF6ElTVl6tqv9nebpJ9klSSy/vad0ry4yS3zsF7rkhyY5KfJHn5gOU/l+QzSR5IcleSv5ntGqS5YHiT5tYLqmob4GnAfwf+bHM3kGTxrFclDTCi37XHJPmFnvkXA/9vjt7rKuD3gcv7FyTZCrgI+CLwBGAPYGJCtLQphjdpBKrqDuB84BcAkjw2yWlJ1iS5I8nbNxzKbA89/d8kpyS5Bzi5/3BUkmcm+XqS77U/n9mz7IlJLm73JlwE7NSzbMPej8Xt/I7tHsHVSe5Ncs50nyHJ7ya5vt3uN5I8rW3fvz3Md197OPConnXOTPLeJJ9t1/tqkie1y05N8rd97/HpJP97pu32rbPRobr2Mz45yXLgJcAb2z2g/9ouvzXJc9vpRyX5u/Y7WN1OP6pdtizJqiRvSLK2Ha9XTFPHrya5pmf+80m+1jN/aZIX9r5/ksOBtwC/3dZ3Vc8m925/Dx5IcmGSnRigp8a3tHuPbk3ykp7lj0ryt0m+3e4BPjXJ1n3rvinJd4AzNrT1rL+p8X1cknOT3N9+1icNqrHPWcDxPfMvAz4wxHqbrareW1VfAH44YPHLgdVV9a6q+s+q+mFVXT0XdUizzfAmjUCSPYEjgCvapvcD64EnA08FDgN6z3F6OnAL8HjgL/q2tSPwWeDvgccB7wI+m+RxbZcPA5fRhLY/5+F/KPudBfwscGD7XqdMU/+LgJNp/tBuBxwF3J1kCvhX4MJ2/T8APpSk97DbccBbgR2Am3s+z4dpQkva99ih/R4+MuR2Z1RVK4APAX/THi5+wYBufwocAhwELAEO5uF7SJ8APBbYHXgV8N621n5fAZ6c5jDgYpqgvkeSbduw9EvAl/vq+xzwl8BH2/qW9Cx+MfCK9vNvBfzRJj7qE2jGe3ea8V7R8139NfDz7ed7ctvnxL51dwT2Bpb3bnSIcXgvTTDaFXhl+5rJB4FjkyxKsj+wLfDVTa2Q5Oo2PA56/cMQ7znIIcCtSc5vQ++Xkvy3LdyWNFKGN2lunZPkPuBS4GLgL5PsAjwfeH37L/61NKHp2J71VlfVu6tqfVU92LfN3wBuqqqz2uVnAzcAL0iyF83h2f9TVT+qqkto/vhuJMmubR0nVNW9VbWuqi6e5nO8miYAfb0aN1fVbTR/ALcB/qqqflxVXwQ+QxPYNvhkVX2tqtbTBKmD2vYvAwU8q50/BvhKVa0ecruz5SXA26pqbVV9lyZovrRn+bp2+bqqOg/4PrBRiKyqHwIrgWcDS4Gracb90Pbz3FRVd29GXWdU1Tfb8f8YP/3eprNhzC+mCfe/1Qbj3wX+V1XdU1UP0ITF3t+1nwAntev2/65NOw5p9hT/T+DE9vf4Wpp/lMxkFXAj8FyaoDnjXreq+sWq2n6a1+8P8Z6D7EHzPfw9sBvNd/bpNIdTpYnmuTTS3HphVX2+t6H91/0UsKbd6QTNP6Ru7+nWO91vN+C2vrbbaPao7AbcW1X/2bdszwHb2RO4p6runelDtH2/NU0tt1fVTwbUssF3eqZ/QBMGqKpK8hGaQHYJzZ6mDeccDbPd2dL/fd7Wtm1wdxs8N/ivzzDAxcAymoByMXAv8CvAj9r5zTHwe5vGoDHfDdiZZs/qZT2/awF6rzb+bhs8B9nUOOxM8zfk9r5lw/gAzWHLZ9KE3X2HXG82PQhcWlXnA7SH8P8M2J/mXDlpYrnnTRq922n+mO/Us/dgu6o6sKdPbWL91TSHuHrtBdwBrAF2SPKYvmXT1bFjku2HrHnQ+UyrgT2T9P6/ZEMtwzgbOCbJ3jSHij+xBdv9T5qAAkCSJ/Qt39R3ueG9er/Pvdq2LbEhvD27nb6YJrz9CtOHt5nqG8agMV8N3EUTUg7s+V17bHsRzTDvv6lx+C7Nof89+5YN4xM0e5BvaffgblJ7rt33p3mdOuR79rua2fnupZEzvEkjVlVraM4hemeS7ZL8TJInJfmVITdxHvDzSV6cZHGS3wYOAD7T/iFcCbw1yVZJfhkYdJ7XhjrOB/4hyQ5JppI8e5r3fB/wR0l+KY0nt4HrqzTh6Y3t+sva9/vIkN/FFTQh4H3ABVV1X7toc7Z7FXBgkoOSPJrm3LxedwI/t4kyzgb+LMnO7UUBJ7LlVx3+O80h1YOBr1XVdTTB8Ok0excHuRPYpy8gbYkNY/4s4Ejg4+0es38GTknyeIAkuyd53pDbnHYc2tvNfJLmgpqfTXIAmz6/8r+0ewmfw8PP89xU/wPbcwIHvU6Ybr32+3g0zd7GqSSP7vmePwgckubCkUXA62nC7vXD1CSNk+FNGo+X0ZyE/g2aQ2v/QnPS94za86aOBN4A3A28ETiyqu5qu7yYJizcA5zEps8peinNOV03AGtp/oANes+P01xo8GHgAeAcYMeq+jHNxQvPp/nD9w/Ay6rqhmE+S+tsmvOfPtzzfkNvt6q+CbwN+DxwE815Zr1OAw5oT24fdDXt22kC79XANTS3ldiiGwC3oeRy4Lr2M0BzIcNt7bmNg3y8/Xl3+u6Bthm+Q/N7tJrmvMITer6rN9FcKPIfSe6n+Z6GuvBjiHF4Hc3h3O8AZwJnDFtwVa2sqkGH4mfThTR7Hp8JrGinn92+/43A7wCn0nx3RwNH9YybNLFS5V5jSeqqdm/YB6tqj3HXImk03PMmSZLUIYY3SZKkDhnZYdMkp9Ocp7O2qjbcZf6j/PTci+2B+6pqo3sZpXnm3QPAQ8D6qlo6kqIlSZImzCjD27Npbm75gQ3hrW/5O4HvVdXbBiy7FVjac0K2JEnSgjSym/RW1SVJ9hm0rL0L+G/RXDouSZKkaUzKExaeBdxZVTdNs7yAC5MU8E/VPK9wRjvttFPts88+s1SiJEnS3LnsssvuqqqdZ+o3KeHtOJp7PU3n0Kpa3d5k8qIkN7TPbNxIkuW0D1fea6+9WLly5exXK0mSNMuSDPWIubFfbZpkMfA/gI9O16d9UDXtTS4/RXP38un6rqiqpVW1dOedZwyvkiRJnTL28EZzZ/UbqmrVoIVJHpNk2w3TwGHAtSOsT5IkaWKMLLwlOZvmMTH7JVmV5FXtomPpO2SaZLck57WzuwCXJrkK+Brw2ar63KjqliRJmiSjvNr0uGnaXz6gbTVwRDt9C7BkTouTJEnqiEk4bCpJkqQhGd4kSZI6xPAmSZLUIYY3SZKkDjG8SZIkdYjhTZIkqUMm5fFYc+KaO77HoX/1Rf74efvxwqfuPrDPOVfcwTsuuJHV9z3Ibttvvcm+m9t/LrdtLQv7c0qSFq5FJ5988rhrmDNvf+e7T66nPJeLv/ld9thha56y63YPW37OFXfwJ5+8hnt+8GMAHvjh+mn7bm7/udy2tSzszylJmp/e+ta3rjn55JNXzNRvQRw2fXDdQ7zjghs3an/HBTfy4LqHhuq7uf3nctvWsrA/pyRpYVsQ4Q1g9X0PDtU2W+1zuW1rWdifU5K0sC2Y8Lbb9lsP1TZb7XO5bWtZ2J9TkrSwLYjwtvXUIv74eftt1P7Hz9uPracWDdV3c/vP5batZWF/TknSwjbvL1h4yrLf5MQXHDDwqr2n7Lode+ywNV+8YS0PVbH79ltP23dz+8/ltq1lYX9OSdL8NOwFC6mqUdQzFjvuvX/dc9v1M/b77X/6CgAffc0zhtru5vSfy21by+i3PWm1SJLmjySXVdXSmfqN7LBpktOTrE1ybU/byUnuSHJl+zpimnUPT3JjkpuTvHlUNUuSJE2aUZ7zdiZw+ID2U6rqoPZ1Xv/CJIuA9wLPBw4AjktywJxWKkmSNKFGFt6q6hLgni1Y9WDg5qq6pap+DHwEOHpWi5MkSeqISbja9HVJrm4Pq+4wYPnuwO0986vaNkmSpAVn3OHtH4EnAQcBa4B3DuiTAW3TXmWRZHmSlUlWrlu3bnaqlCRJmhBjDW9VdWdVPVRVPwH+meYQab9VwJ4983sAqzexzRVVtbSqlk5NTc1uwZIkSWM21vCWZNee2d8Erh3Q7evAvkmemGQr4Fjg3FHUJ0mSNGkWj+qNkpwNLAN2SrIKOAlYluQgmsOgtwKvafvuBryvqo6oqvVJXgdcACwCTq+q60ZVtyRJ0iQZWXirquMGNJ82Td/VwBE98+cBG91GRJIkaaEZ9wULkiRJ2gyGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqEMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqkJGFtySnJ1mb5NqetnckuSHJ1Uk+lWT7ada9Nck1Sa5MsnJUNUuSJE2aUe55OxM4vK/tIuAXquoXgW8Cf7KJ9X+1qg6qqqVzVJ8kSdLEG1l4q6pLgHv62i6sqvXt7H8Ae4yqHkmSpC6apHPeXgmcP82yAi5MclmS5SOsSZIkaaIsHncBAEn+FFgPfGiaLodW1eokjwcuSnJDuydv0LaWA8sBttn1SXNSryRJ0riMfc9bkuOBI4GXVFUN6lNVq9ufa4FPAQdPt72qWlFVS6tq6dTU1FyULEmSNDZjDW9JDgfeBBxVVT+Yps9jkmy7YRo4DLh2UF9JkqT5bpS3Cjkb+AqwX5JVSV4FvAfYluZQ6JVJTm377pbkvHbVXYBLk1wFfA34bFV9blR1S5IkTZKRnfNWVccNaD5tmr6rgSPa6VuAJXNYmiRJUmeM/Zw3SZIkDc/wJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqEMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqkJGGtySnJ1mb5Nqeth2TXJTkpvbnDtOse3zb56Ykx4+uakmSpMkx6j1vZwKH97W9GfhCVe0LfKGdf5gkOwInAU8HDgZOmi7kSZIkzWcjDW9VdQlwT1/z0cD72+n3Ay8csOrzgIuq6p6quhe4iI1DoCRJ0rw3Cee87VJVawDan48f0Gd34Pae+VVt20aSLE+yMsnKdevWzXqxkiRJ4zQJ4W0YGdBWgzpW1YqqWlpVS6empua4LEmSpNGahPB2Z5JdAdqfawf0WQXs2TO/B7B6BLVJkiRNlBnDW3s16Eyv7R9BDecCG64ePR749IA+FwCHJdmhvVDhsLZNkiRpQVk8RJ/V7WvQocsNFgF7zbShJGcDy4CdkqyiuYL0r4CPJXkV8G3gRW3fpcAJVfXqqronyZ8DX2839baq6r/wQZIkad4bJrxdX1VP3VSHJFcM82ZVddw0i35tQN+VwKt75k8HTh/mfSRJkuarYc55ewZAkrf3L0iyqLePJEmS5taM4a2qfthO7p7kxRvakzwe+HxfH0mSJM2hYQ6bbvAa4IIkN9PcpuMM4E1zUpUkSZIGmjG8JfkAcDlwBfBa4MPAeuCFVXXz3JYnSZKkXsOc8/b+tt8raYLbPsC9wO8kOWbuSpMkSVK/Gfe8VdUXaB4YD0CSxcABwBLgEOBf5qw6SZIkPczmnPMGQFWtB65uX2fNekWSJEma1jBPWLh8NvpIkiTpkRtmz9v+Sa7exPIAj52leiRJkrQJw4S3pwzR56FHWogkSZJmNswFC7cBJPk88IaqumrOq5IkSdJAw9wqZIM3AqckOSPJrnNVkCRJkqY3dHirqsur6jnAZ4DPJTkpydZzV5okSZL6bc6eN5IEuBH4R+APgJuSvHQuCpMkSdLGhg5vSS4F7gBOAXYHXg4sAw5OsmJLC0iyX5Ire173J3l9X59lSb7X0+fELX0/SZKkLtucm/SeAFxXVdXX/gdJrt/SAqrqRuAggCSLaALipwZ0/XJVHbml7yNJkjQfDB3equraTSz+jVmoBeDXgG9tuMJVkiRJD7dZ57xNp6pumY3tAMcCZ0+z7BlJrkpyfpIDp9tAkuVJViZZuW7dulkqS5IkaTLMSnibDUm2Ao4CPj5g8eXA3lW1BHg3cM5026mqFVW1tKqWTk1NzU2xkiRJYzIx4Q14PnB5Vd3Zv6Cq7q+q77fT5wFTSXYadYGSJEnjNknh7TimOWSa5AntbUpIcjBN3XePsDZJkqSJsDlXm86ZJD8L/Drwmp62EwCq6lTgGOD3kqwHHgSOHXDVqyRJ0rw3EeGtqn4APK6v7dSe6fcA7xl1XZIkSZNmkg6bSpIkaQaGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqEMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqkIkJb0luTXJNkiuTrBywPEn+PsnNSa5O8rRx1ClJkjROi8ddQJ9fraq7pln2fGDf9vV04B/bn5IkSQvGxOx5G8LRwAeq8R/A9kl2HXdRkiRJozRJ4a2AC5NclmT5gOW7A7f3zK9q2x4myfIkK5OsXLdu3RyVKkmSNB6TFN4Oraqn0RwefW2SZ/ctz4B1aqOGqhVVtbSqlk5NTc1FnZIkSWMzMeGtqla3P9cCnwIO7uuyCtizZ34PYPVoqpMkSZoMExHekjwmybYbpoHDgGv7up0LvKy96vQQ4HtVtWbEpUqSJI3VpFxtugvwqSTQ1PThqvpckhMAqupU4DzgCOBm4AfAK8ZUqyRJ0thMRHirqluAJQPaT+2ZLuC1o6xLkiRp0kzEYVNJkiQNx/AmSZLUIYY3SZKkDjG8SZIkdYjhTZIkqUMMb5IkSR1ieJMkSeoQw5skSVKHGN4kSZI6xPAmSZLUIYY3SZKkDjG8SZIkdYjhTZIkqUMMb5IkSR0y9vCWZM8k/5bk+iTXJfnDAX2WJflekivb14njqFWSJGncFo+7AGA98IaqujzJtsBlSS6qqm/09ftyVR05hvokSZImxtj3vFXVmqq6vJ1+ALge2H28VUmSJE2msYe3Xkn2AZ4KfHXA4mckuSrJ+UkOHGlhkiRJE2JiwluSbYBPAK+vqvv7Fl8O7F1VS4B3A+dsYjvLk6xMsnLdunVzV7AkSdIYTER4SzJFE9w+VFWf7F9eVfdX1ffb6fOAqSQ7DdpWVa2oqqVVtXRqampO65YkSRq1sYe3JAFOA66vqndN0+cJbT+SHExT992jq1KSJGkyTMLVpocCLwWuSXJl2/YWYC+AqjoVOAb4vSTrgQeBY6uqxlGsJEnSOI09vFXVpUBm6PMe4D2jqUiSJGlyjf2wqSRJkoZneJMkSeoQw5skSVKHGN4kSZI6xPAmSZLUIYY3SZKkDjG8SZIkdYjhTZIkqUMMb5IkSR1ieJMkSeoQw5skSVKHGN4kSZI6xPAmSZLUIYY3SZKkDpmI8Jbk8CQ3Jrk5yZsHLH9Uko+2y7+aZJ/RVylJkjR+Yw9vSRYB7wWeDxwAHJfkgL5urwLuraonA6cAfz3aKiVJkibD2MMbcDBwc1XdUlU/Bj4CHN3X52jg/e30vwC/liQjrFGSJGkipKrGW0ByDHB4Vb26nX8p8PSqel1Pn2vbPqva+W+1fe7a1Lb33W77+vzRR81YwzfW3A/AAbtuN1TNm9N/LrdtLaPf9ihq+c7Oe/KKs989VH9J0vyR5LKqWjpTv8WjKGYGg/ag9SfKYfo0HZPlwPJ29kf7fPCsax9BbZosOwGbDOzzxSs/8p5xlzAKC2Y8FwjHc/5xTEdv72E6TUJ4WwXs2TO/B7B6mj6rkiwGHgvcM2hjVbUCWAGQZOUwCVbd4HjOL47n/OJ4zj+O6eSahHPevg7sm+SJSbYCjgXO7etzLnB8O30M8MUa9/FeSZKkMRj7nreqWp/kdcAFwCLg9Kq6LsnbgJVVdS5wGnBWkptp9rgdO76KJUmSxmfs4Q2gqs4DzutrO7Fn+ofAi7Zg0yseYWmaLI7n/OJ4zi+O5/zjmE6osV9tKkmSpOFNwjlvkiRJGtK8DG8zPW5Lky/J6UnWtvf429C2Y5KLktzU/txhnDVqeEn2TPJvSa5Pcl2SP2zbHdMOSvLoJF9LclU7nm9t25/YPsLwpvaRhluNu1YNL8miJFck+Uw773hOqHkX3oZ83JYm35nA4X1tbwa+UFX7Al9o59UN64E3VNX+wCHAa9v/Lh3TbvoR8JyqWgIcBBye5BCaRxee0o7nvTSPNlR3/CFwfc+84zmh5l14Y7jHbWnCVdUlbHwvv97HpL0feOFIi9IWq6o1VXV5O/0AzR+I3XFMO6ka329np9pXAc+heYQhOJ6dkmQP4DeA97XzwfGcWPMxvO0O3N4zv6ptU/ftUlVroAkDwOPHXI+2QJJ9gKcCX8Ux7az2ENuVwFrgIuBbwH1Vtb7t4v97u+XvgDcCP2nnH4fjObHmY3gb+lFakkYryTbAJ4DXV9X9465HW66qHqqqg2ieinMwsP+gbqOtSlsiyZHA2qq6rLd5QFfHc0JMxH3eZtkwj9tSN92ZZNeqWpNkV5p/8asjkkzRBLcPVdUn22bHtOOq6r4kX6I5l3H7JIvbvTX+v7c7DgWOSnIE8GhgO5o9cY7nhJqPe96GedyWuqn3MWnHA58eYy3aDO35M6cB11fVu3oWOaYdlGTnJNu301sDz6U5j/HfaB5hCI5nZ1TVn1TVHlW1D83fzC9W1UtwPCfWvLxJb/uvh7/jp4/b+osxl6TNlORsYBmwE3AncBJwDvAxYC/g28CLqqr/ogZNoCS/DHwZuIafnlPzFprz3hzTjknyizQnsC+i2Qnwsap6W5Kfo7lIbEfgCuB3qupH46tUmyvJMuCPqupIx3NyzcvwJkmSNF/Nx8OmkiRJ85bhTZIkqUMMb5IkSR1ieJMkSeoQw5skSVKHGN4kSZI6xPAmSZLUIYY3SRpSkucmOWvcdUha2AxvkjS8JTR3mpeksTG8SdLwlgBXJHlUkjOT/GX73FZJGpnF4y5AkjpkCbAWuAB4X1V9cMz1SFqAfLapJA0hyRRwF3Ab8Jqq+sqYS5K0QHnYVJKGcwDwdWA98NCYa5G0gBneJGk4S4B/B44Fzkiyy5jrkbRAGd4kaThLgGur6pvAm4CPtYdSJWmkPOdNkiSpQ9zzJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqEMObJElShxjeJEmSOuT/A1C1vl1FwnYbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAADiCAYAAABwdKKfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHChJREFUeJzt3XuYXXV56PHvyxBkACECUcmABq1GOAWMJ0dBtKVKnyBQpR56REXFS8GnwrGtBom1Xmq52FhRjz3WiOAFiWDMyfGoJV6Qo1QLBgYJyhmL3CdcgjJAdZQQ3vPHWkN2xj0ze09mr337fp5nnuy9Lr/17nX7vVnrt34rMhNJkiRVY6d2ByBJktRPTL4kSZIqZPIlSZJUIZMvSZKkCpl8SZIkVcjkS5IkqUImX9IkEfHiiBhpdxydKCKujIi37MD8/xERz5jLmMpyb4uIo+e63B21I/tSRJwSEVfVfG/JuqtSRLw/Ii7egfn/JSLeMJcxSe1g8qW+NVWFnZnfz8zF7Yipl9RL1DJzj8y8pV0xVW0u96V+W3f1ErXMfFlmfq5dMUlzxeRL6hARsXO7Y9DccXtKmorJlzRJRBwVEXfVfL8tIt4ZETdExIMRcWlE7Foz/viIuD4ixiLiBxFxaM24syLi5xHxcET8NCL+tGbcKRHxrxFxfkT8Enh/nVgGIuLdNWVcGxEHlONeGBE/KmP6UUS8sGa+KyPig2X5D0fENyNi33Lc5RFx+qTl/DgiXjlTuZPm2e7KREQsioiMiJ0j4mzgxcAnyttlnyinyYj4vfLzXhHx+YjYHBG3R8R7ImKnmnVzVUR8OCIeiIhbI+JlM2y6/1Ku4wci4qKJbRQRN0bEn9TEOS8i7o+I59b5TUdFxF3lOr+/3PavrRn/hDKmOyLi3oj454gYnDTvuyLiHuCiOvvSQeW2GYuIn0TEy2vG7RMRX42IhyLiGuCZk2KrXXeDEfGP5Xp7sFxXg1Nsp1eU++dD5X50TDl8Ybm8X0bEzRHx55O27WXl9nm4jHVpOe6siFgzaRkfi4iPz1RuvXU9adhtEXF0GeO7gVeV+8+Py/GPX02NiJ3Kfeb2iLivjHWvctzEvviGclvdHxF/Uy8OqR1MvqTG/DfgGOBA4FDgFICIeB5wIXAasA/wKeCrEfGEcr6fUyQhewEfAC6OiP1qyn0BcAvwZODsOsv9a+DVwLHAnsCbgF9HxN7A14GPl8v9CPD1iNinZt7XAG8sy94FeGc5/JKyTMrfcDDw9HL+RsqdUWb+DfB94PTydtnpdSb7HxTr5RnAHwKvL+Od8AJgBNgX+AfgMxER0yz2tcAyiqTl2cB7yuGfB06ume5Y4O7MvH6Kcp5aLnMIeAOwKiImbh1+qCz7ucDvldO8d9K8e1Osz1NrC42IecD/Ab5JsU3OAL5YU/Y/Ab8B9qPYzm+a5rd+GPjPwAvL5Z0JPDZ5ooh4fvn7lwPzgT8AbitHrwbuAhYCJwLnRMRLa2Z/OfClcr6vAp+ome/YiNizXMYAxfFxSYPlzigzLwfOAS4t95/D6kx2Svn3RxT70B41MU54EbAYeCnw3og4qJk4pJbJTP/868s/ikro6DrDjwLumjTdyTXf/wH45/LzJ4EPTpp/BPjDKZZ5PfCK8vMpwB0zxDgyMf2k4a8Drpk07IfAKeXnK4H31Iz7C+Dy8vMTgV8BTy+/nw1c2ES5byk/vx+4uGa6RUACO0+etmaapEhaBoDfAgfXjDsNuLJm3dxcM263ct6nTrMt31rz/Vjg5+XnhcDDwJ7l9zXAmVOUcxTwKLB7zbDLgL8Folxvz6wZdwRwa828jwC71tuXKJLwe4CdasavLtfjALAFeE7NuHOAq+qsu52AceCwBvbxTwHn1xl+ALAVeGLNsHOBz9Zs22/XjDsYGK/5fhXw+vLzH9es60bKvbjecTb5mGTS/lVn//sO8Bc14xaX63Bntu2L+9eMvwY4aaZ15p9/Vfx55UtqzD01n39N8b9sKK5wvKO8jTQWEWMUFdBCgIh4fWy7JTkG/D7FVZUJd86w3AMorp5NthC4fdKw2ymuxEwbc2Y+THF166Ry3EnAF5sody7sS3E1rnZZU8afmb8uP+7B1GrX5e2U2yAzNwH/CvzXiJgPvIxtv7eeBzLzV3XKWkCRBF5bsz0vL4dP2JyZv5mi3IXAnZlZe4Vq4jcvoEgaJv+GevYFdqX+fjHZdPvPL8t9YXIsEybvP7vGtnZstVdPX8O2q16NlDtXJu+rt1Osw6fUDJvquJXayuRL2jF3Amdn5vyav90yc3VEPB34NHA6sE9mzgdupLiCMiEbKP+ZdYZvokj8aj0NGG0w7tXAqyPiCGAQ+O4syv0VRTIy4amTxk/32+6nuEpRu6xm4q/ngEllbar5/jmKW49/BvwwM6dbzpMiYvc6Zd1PccXpP9Vs670ys7ZCn+43bwIOiLJdW03Zo8Bmiituk39DPfdT3J6st19MNt3+s3dEPLFOLI34MnBUROwP/Cnbkq9myt1u/ylvX9YmsjMdG5P31adRrMN7G/kBUjuZfKnfzYuIXWv+mn1C7dPAWyPiBVHYPSKOKyuf3SkqkM0AEfFGiitfzbgA+GBEPKss/9Cy/dU3gGdHxGuiaOD+KopbQ19rsNxvUFRcf0fRruaxmuGNlns98AcR8bSyofOKSePvpWiL8zsycyvF7byzI+KJZaL618Cs+4AC3hYR+5ft1t4NXFozbh3wPODtFG2gZvKBiNglIl4MHA98uVxHnwbOj4gnA0TEUEQsazC+qykSjjOjaPR/FPAnwJfK9bEWeH9E7Fa2w6vbn1UZx4XAR8rG7QMRcURNO8NanwHeGBEvLRuoD0XEczLzTuAHwLnlfn8o8GamvyJYG8NmiluAF1Hcdr2pHN5MuT+juJp2XNke7j1A7W+4F1g0KVmttRr4q4g4MCL2YFsbsUcb+Q1SO5l8qd99g+JqxsTf+5uZOTM3AH9O0dD3AeBmysb4mflT4B8p2kzdCxxCcfurGR+hSFK+CTxEUZkOZuYvKJKCdwC/oGhwfXxm3t9g3L+lqOyPZttVC5opNzO/RZHg3ABcy+8maB8DTozi6cOP1wnjDIpk5BaKNkSXUCQVs3UJxXq6pfz7+5pYx4GvUDwwsXaGcu6h2JabKJKGt2bm/yvHvYtiG/9bRDwEfJuirdGMMvMRikbsL6O4evU/KdpNTZR9OsVtsXuAz1IkNlN5J7AR+BHwS4oHAX7nfJ6Z11A8xHA+8CDwf9l2tejVFG2jNgH/C3hfuU0bdQmT9p9mys3MBynaIl5AcWXsVxQN9Sd8ufz3FxFxXZ3lXwh8AfgecCvF1cAzmohfapvInOnKriR1v4h4L/DszDx5mmmOomjkvX9lgUnqO3YCKKnnlbci30zxNKcktZW3HSX1tCg6+bwT+JfM/F6745EkbztKkiRVyCtfkiRJFTL5kiRJqlBHN7jfd999c9GiRe0OQ5IkaUbXXnvt/Zm5YKbpOjr5WrRoERs2bGh3GJIkSTOKiKleC7YdbztKkiRVyORLkiSpQiZfkiRJFTL5kiRJqpDJlyRJUoVMviRJkipk8iVJklQhky9JkqQKmXxJkiRVqNLkKyL+KiJ+EhE3RsTqiNi1yuVLkiS1W2WvF4qIIeC/Awdn5nhEXAacBHy2qhikVlo3PMrK9SNsGhtn4fxBli9bzAlLhuZk+m4t21j6+3dKqq/qdzvuDAxGxBZgN2BTxcuXGtZsBbZi7UbGt2wFYHRsnBVrNwLUnaeZ6bu1bGPp799ZO4/JmrS9ym47ZuYo8GHgDuBu4MHM/GZVy5eaMVHJjI6Nk2yrZNYNj9adfuX6kccrpAnjW7aycv3IDk/frWUbS3//Tmj+OJL6RWXJV0Q8CXgFcCCwENg9Ik6uM92pEbEhIjZs3ry5qvCk7TRbyWwaG2/Z8G4t21j6+3dC88eR1C+qbHB/NHBrZm7OzC3AWuCFkyfKzFWZuTQzly5YsKDC8NQP1g2PcuR5V3DgWV/nyPOumPJ/4M1WMgvnD7ZseLeWbSz9/Tuh+eOo0eNT6nZVJl93AIdHxG4REcBLgZsqXL76XDO3QJqtZJYvW8zgvIHthg3OG2D5ssU7PH23lm0s/f07obnjyFuU6idVtvm6GlgDXAdsLJe9qqrlS83cAmm2kjlhyRDnvvIQdhkoDqmh+YOc+8pDpmxY3Mz03Vq2sfT374TmjiNvUaqfRGa2O4YpLV26NDds2NDuMNQjDjzr69Tb2wO49bzjfmf4uuFRzlxzA49sfYyhBp/SetWnfgjApacd0VBMzUzfrWUbS/Vld1IsjR5HzR6fUieKiGszc+lM01Xd1YTUNgvnDzJap63JVLdGTlgyxOpr7gAar5Qkba/R46jZ41PqZr5eSH2j2VuJkqrj8al+4pUvdbVmOnCcGN7srURJrTeb49MOXNWtTL7UtWbT27a3EqXO1czxOZvjX+oU3nZU1/LpKKl/efyrm5l8qWs124GjpN7h8a9uZvKlrtVsR6iSeofHv7qZyZe6lk9HSf3L41/dzAb36lo+vSj1L49/dTOTL3U1n16U+pfHv7qVyZc6iv32SGoFzy3qJCZf6hj22yOpFTy3qNPY4F4dw357JLWC5xZ1GpMvdQz77ZHUCp5b1GlMvtQx7LdHUit4blGnMflSx7DfHkmt4LlFncYG9+oY9tsjqRU8t6jTmHypo9hvj6RW8NyiTuJtR0mSpAqZfEmSJFXI5EuSJKlCtvlSy/laD0ndxvOWWsnkSy3laz0kdRvPW2o1bzuqpXyth6Ru43lLrWbypZbytR6Suo3nLbWayZdaytd6SOo2nrfUaiZfailf6yGp23jeUqvZ4F4t5Ws9JHUbz1tqtUqTr4iYD1wA/D6QwJsy84dVxqDq+VoPSd3G85ZaqeorXx8DLs/MEyNiF2C3ipcvSZLUVpUlXxGxJ/AHwCkAmfkI8EhVy5ckSeoEVTa4fwawGbgoIoYj4oKI2L3C5UuSJLVdlcnXzsDzgE9m5hLgV8BZkyeKiFMjYkNEbNi8eXOF4UmSJLVelW2+7gLuysyry+9rqJN8ZeYqYBXA0qVLs7rw1CjfeSZJ2/O8qGZUlnxl5j0RcWdELM7MEeClwE+rWr7mhu88k6TteV5Us6ruZPUM4IsRcQPwXOCcipevHeQ7zyRpe54X1axKu5rIzOuBpVUuU3PLd55J0vY8L6pZvl5ITfGdZ5K0Pc+LapbJl5riO88kaXueF9Us3+2opvjOM0nanudFNcvkS03znWeStD3Pi2qGtx0lSZIqZPIlSZJUIZMvSZKkCpl8SZIkVcjkS5IkqUImX5IkSRUy+ZIkSaqQ/XwJgHXDo6xcP8KmsXEW2kGgJLWM51uZfIl1w6OsWLuR8S1bARgdG2fF2o0AnhAkaQ55vhV421HAyvUjj58IJoxv2crK9SNtikiSepPnW4HJl4BNY+NNDZckzY7nW4HJl4CF8webGi5Jmh3Pt4IGkq+I2LuBv/lVBKvWWL5sMYPzBrYbNjhvgOXLFrcpIknqTZ5vBY01uN9U/sU00wwAT5uTiFS5iUaeZ665gUe2PsaQT99IUkt4vhU0lnzdlJlLppsgIobnKB61yQlLhlh9zR0AXHraEW2ORpJ6l+dbNdLm6wiAiPj7ySMiYqB2GkmSJE1vxuQrM39TfhyKiNdMDI+IJwPfnjSNJEmSptFMJ6unAesj4mYggYuAd7UkKkmSpB41Y/IVEZ8HrgOGgbcBlwCPAidk5s2tDU+SJKm3NNLm63PldG+iSLwWAQ8AJ0fEia0LTZIkqffMeOUrM78DfGfie0TsDBwMHAYcDqxpWXSSJEk9pukXa2fmo8AN5d8X5jwiSZKkHtZID/fXzcU0kiRJauzK10ERccM04wPYa47ikSRJ6mmNJF/PaWCarY0usOyYdQMwmpnHNzqfJElSL2ikwf3tABHxbeAdmfnjHVzm24GbgD13sBxNY93wKCvXj7BpbJyFvjtMkrqW5/Pe00hXExPOBM6PiIsiYr/ZLCwi9geOAy6YzfxqzLrhUVas3cjo2DgJjI6Ns2LtRtYNj7Y7NElSEzyf96aGk6/MvC4zXwJ8Dbg8It4XEYNNLu+jFEncY03OpyasXD/C+Jbt7wSPb9nKyvUjbYpIkjQbns97UzNXvoiIAEaATwJnAP8eEa9rcN7jgfsy89oZpjs1IjZExIbNmzc3E55Km8bGmxouSepMns97U8PJV0RcBYwC5wNDwCnAUcDzI2JVA0UcCbw8Im4DvgS8JCIunjxRZq7KzKWZuXTBggWNhqcaC+fXvyA51XBJUmfyfN6bmrny9VZgKDP/ODP/NjO/lpk3Z+YZwItnmjkzV2Tm/pm5CDgJuCIzT55d2JrO8mWLGZw3sN2wwXkDLF+2uE0RSZJmw/N5b2q4h/vMvHGa0cfNQSyaIxNPwZy55gYe2foYQz4dI0ldyfN5b2r69UL1ZOYtTU5/JXDlXCxb9Z2wZIjV19wBwKWnHdHmaCRJs+X5vPc01eBekiRJO8bkS5IkqUImX5IkSRUy+ZIkSaqQyZckSVKFTL4kSZIqZPIlSZJUIZMvSZKkCpl8SZIkVWhOerhX660bHmXl+hE2jY2z0NdLSJKmYH3R+Uy+usC64VFWrN3I+JatAIyOjbNi7UYADyhJ0uOsL7qDtx27wMr1I48fSBPGt2xl5fqRNkUkSepE1hfdweSrC2waG29quCSpP1lfdAeTry6wcP5gU8MlSf3J+qI7mHx1geXLFjM4b2C7YYPzBli+bHGbIpIkdSLri+5gg/suMNFI8sw1N/DI1scY8ukVSVId1hfdweSrS5ywZIjV19wBwKWnHdHmaCRJncr6ovN521GSJKlCJl+SJEkVMvmSJEmqkMmXJElShUy+JEmSKmTyJUmSVCGTL0mSpArZz1cbrRseZeX6ETaNjbPQjvAkSRWzHmoPk682WTc8yoq1Gx9/+/zo2Dgr1m4EcMeXJLWc9VD7eNuxTVauH3l8h58wvmUrK9ePtCkiSVI/sR5qH5OvNtk0Nt7UcEmS5pL1UPtUlnxFxAER8d2IuCkifhIRb69q2Z1o4fzBpoZLkjSXrIfap8orX48C78jMg4DDgbdFxMEVLr+jLF+2mMF5A9sNG5w3wPJli9sUkSSpn1gPtU9lDe4z827g7vLzwxFxEzAE/LSqGDrJRGPGM9fcwCNbH2PIp0wkSRWyHmqftjztGBGLgCXA1e1Yfqc4YckQq6+5A4BLTzuizdFIkvqN9VB7VN7gPiL2AL4C/GVmPlRn/KkRsSEiNmzevLnq8CRJklqq0uQrIuZRJF5fzMy19abJzFWZuTQzly5YsKDK8CRJklquyqcdA/gMcFNmfqSq5UqSJHWSKq98HQm8DnhJRFxf/h1b4fIlSZLarsqnHa8CoqrlSZIkdSLf7TiHfEGpJKmXWc/NDZOvOeILSiVJvcx6bu74bsc54gtKJUm9zHpu7ph8zRFfUCpJ6mXWc3PH5GuO+IJSSVIvs56bOyZfc8QXlEqSepn13Nyxwf0c8QWlkqReZj03d0y+5pAvKJUk9TLrubnhbUdJkqQKmXxJkiRVyNuO07AnX0mSZsc6dGomX1OwJ19JkmbHOnR63nacgj35SpI0O9ah0zP5moI9+UqSNDvWodMz+ZqCPflKkjQ71qHTM/magj35SpI0O9ah07PB/RTsyVeSpNmxDp2eydc07MlXkqTZsQ6dWt8lX/Y7IklS5+mn+rmvki/7HZEkqfP0W/3cVw3u7XdEkqTO02/1c18lX/Y7IklS5+m3+rmvki/7HZEkqfP0W/3cV8mX/Y5IktR5+q1+7voG9808HWG/I5IkdZ7Z1M/d/HRkVydfs3k6wn5HJEnqPM3Uz93+dGRX33bst6cjJElS99f/XZ189dvTEZIkqfvr/0qTr4g4JiJGIuLmiDhrpuk3jj7Ikeddwbrh0brj++3pCEmS1Hz9v254lCPPu4IDz/r6tHlFVSpLviJiAPgn4GXAwcCrI+LgmeabuI9bb0X129MRkiSpufp/on3Y6Ng4yfR5RVWqvPL1fODmzLwlMx8BvgS8opEZp7qPe8KSIc595SHsMlD8jKH5g5z7ykO6orGdJEmanWbq/05sHxaZWc2CIk4EjsnMt5TfXwe8IDNPn2qeZ+6xV55z2Ise/374M/apO91P734IgIP327OhWJqZvpVlG0v1ZRtL9WUbS/VlG0v1ZRtL9WU3Ov2/3fKLxz/fstcQnzq0uO4TwK3nHdfQchoVEddm5tKZpquyq4moM+x3Mr+IOBU4FWCnwT15/S0/Kybc+ugjW35w28aWRqhW2xe4v91BaM64PXuL27P3uE2BeQsWHRIDO+9SfPsZDH8XKPKK+NDxc51XPL2RiapMvu4CDqj5vj+wafJEmbkKWAUQERt+++sHZ8wg1R0iYkMj/yNQd3B79ha3Z+9xm3auKtt8/Qh4VkQcGBG7ACcBX61w+ZIkSW1X2ZWvzHw0Ik4H1gMDwIWZ+ZOqli9JktQJKn29UGZ+A/hGE7OsalUsagu3Z29xe/YWt2fvcZt2qMqedpQkSVKXv15IkiSp23Rk8tXsa4jUeSLiwoi4LyJurBm2d0R8KyL+vfz3Se2MUY2LiAMi4rsRcVNE/CQi3l4Od5t2oYjYNSKuiYgfl9vzA+XwAyPi6nJ7Xlo+HKUuEREDETEcEV8rv7s9O1THJV+zfQ2ROs5ngWMmDTsL+E5mPgv4Tvld3eFR4B2ZeRBwOPC28rh0m3an3wIvyczDgOcCx0TE4cCHgPPL7fkA8OY2xqjmvR24qea727NDdVzyxQ68hkidIzO/B/xy0uBXAJ8rP38OOKHSoDRrmXl3Zl5Xfn6Y4gQ/hNu0K2XhP8qv88q/BF4CrCmHuz27SETsDxwHXFB+D9yeHasTk68h4M6a73eVw9T9npKZd0NRmQNPbnM8moWIWAQsAa7Gbdq1yltU1wP3Ad8Cfg6MZeaj5SSee7vLR4EzgcfK7/vg9uxYnZh8NfQaIknVi4g9gK8Af5mZD7U7Hs1eZm7NzOdSvG3k+cBB9SarNirNRkQcD9yXmdfWDq4zqduzQ1Taz1eDGnoNkbrSvRGxX2beHRH7UfyPW10iIuZRJF5fzMy15WC3aZfLzLGIuJKiLd/8iNi5vFriubd7HAm8PCKOBXYF9qS4Eub27FCdeOXL1xD1rq8Cbyg/vwH4322MRU0o2498BrgpMz9SM8pt2oUiYkFEzC8/DwJHU7Tj+y5wYjmZ27NLZOaKzNw/MxdR1JlXZOZrcXt2rI7sZLXM3j/KttcQnd3mkNSkiFgNHAXsC9wLvA9YB1wGPA24A/izzJzcKF8dKCJeBHwf2Mi2NiXvpmj35TbtMhFxKEUD7AGK/4Rflpl/FxHPoHjIaW9gGDg5M3/bvkjVrIg4CnhnZh7v9uxcHZl8SZIk9apOvO0oSZLUs0y+JEmSKmTyJUmSVCGTL0mSpAqZfEmSJFXI5EuSJKlCJl+SJEkVMvmS1Dci4uiI+EK745DU30y+JPWTwyh6+paktjH5ktRPDgOGI+IJEfHZiDinfG+lJFVm53YHIEkVOgy4D1gPXJCZF7c5Hkl9yHc7SuoLETEPuB+4HTgtM3/Y5pAk9SlvO0rqFwcDPwIeBba2ORZJfczkS1K/OAz4AXAScFFEPKXN8UjqUyZfkvrFYcCNmfkz4F3AZeWtSEmqlG2+JEmSKuSVL0mSpAqZfEmSJFXI5EuSJKlCJl+SJEkVMvmSJEmqkMmXJElShUy+JEmSKmTyJUmSVKH/D9pZrA7Qo3KrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal as sig\n",
    "\n",
    "L = 32  # length of signal x[k]\n",
    "N = 16  # length of signal h[k]\n",
    "M = 16  # periodicity of periodic convolution\n",
    "\n",
    "\n",
    "def periodic_summation(x, N):\n",
    "    \"Zero-padding to length N or periodic summation with period N.\"\n",
    "    M = len(x)\n",
    "    rows = int(np.ceil(M/N))\n",
    "    \n",
    "    if (M < int(N*rows)):\n",
    "        x = np.pad(x, (0, int(N*rows-M)), 'constant')\n",
    "    \n",
    "    x = np.reshape(x, (rows, N))\n",
    "    \n",
    "    return np.sum(x, axis=0)\n",
    "\n",
    "\n",
    "def periodic_convolve(x, y, P):\n",
    "    \"Periodic convolution of two signals x and y with period P.\"\n",
    "    x = periodic_summation(x, P)\n",
    "    h = periodic_summation(y, P)\n",
    "\n",
    "    return np.array([np.dot(np.roll(x[::-1], k+1), h) for k in range(P)], float)\n",
    "\n",
    "\n",
    "# generate signals\n",
    "x = np.ones(L)\n",
    "h = sig.triang(N)\n",
    "\n",
    "# linear convolution\n",
    "y1 = np.convolve(x, h, 'full')\n",
    "# periodic convolution\n",
    "y2 = periodic_convolve(x, h, M)\n",
    "# linear convolution via periodic convolution\n",
    "xp = np.append(x, np.zeros(N-1))\n",
    "hp = np.append(h, np.zeros(L-1))\n",
    "y3 = periodic_convolve(xp, hp, L+N-1)\n",
    "\n",
    "# plot results\n",
    "def plot_signal(x):\n",
    "    plt.figure(figsize = (10, 3))\n",
    "    plt.stem(x)\n",
    "    plt.xlabel(r'$k$')\n",
    "    plt.ylabel(r'$y[k]$')\n",
    "    plt.axis([0, N+L, 0, 1.1*x.max()])\n",
    "\n",
    "\n",
    "plot_signal(x)\n",
    "plt.title('Signal $x[k]$')\n",
    "plot_signal(y1)\n",
    "plt.title('Linear convolution')\n",
    "plot_signal(y2)\n",
    "plt.title('Periodic convolution with period M = %d' %M)\n",
    "plot_signal(y3)\n",
    "plt.title('Linear convolution by periodic convolution');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Change the lengths `L`, `N` and `M` and check how the results for the different convolutions change"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Fast Convolution\n",
    "\n",
    "Using the above derived equality of the linear and periodic convolution one can express the linear convolution $y[k] = x_L[k] * h_N[k]$ by the DFT as\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = \\text{IDFT}_M \\{ \\; \\text{DFT}_M\\{ x_M[k] \\} \\cdot \\text{DFT}_M\\{ h_M[k] \\} \\; \\}\n",
    "\\end{equation}\n",
    "\n",
    "This operation requires three DFTs of length $M$ and $M$ complex multiplications. On first sight this does not seem to be an improvement, since one DFT/IDFT requires $M^2$ complex multiplications and $M \\cdot (M-1)$ complex additions. The overall numerical complexity is hence in the order of $\\mathcal{O}(M^2)$. The DFT can be realized efficiently by the [fast Fourier transformation](https://en.wikipedia.org/wiki/Fast_Fourier_transform) (FFT), which lowers the computational complexity to $\\mathcal{O}(M \\log_2 M)$. The resulting algorithm is known as *fast convolution* due to its computational efficiency. \n",
    "\n",
    "The fast convolution algorithm is composed of the following steps\n",
    "\n",
    "1. Zero-padding of the two input signals $x_L[k]$ and $h_N[k]$ to at least a total length of $M \\geq N+L-1$\n",
    "\n",
    "2. Computation of the DFTs $X[\\mu]$ and $H[\\mu]$ using a FFT of length $M$\n",
    "\n",
    "3. Multiplication of the spectra $Y[\\mu] = X[\\mu] \\cdot H[\\mu]$\n",
    "\n",
    "4. Inverse DFT of $Y[\\mu]$ using an inverse FFT of length $M$\n",
    "\n",
    "The overall complexity depends on the particular implementation of the FFT. Many FFTs are most efficient for lengths which are a power of two. It therefore can make sense, in terms of computational complexity, to choose $M$ as a power of two instead of the shortest possible length $N+L-1$. For real valued signals $x[k] \\in \\mathbb{R}$ and $h[k] \\in \\mathbb{R}$ the computational complexity can be reduced significantly by using a real valued FFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Fast convolution\n",
    "\n",
    "The implementation of the fast convolution algorithm is straightforward. Most implementations of the FFT include zero-padding to a given length $M$, e.g in `numpy` by `numpy.fft.fft(x, M)`. In the following example an implementation of the fast convolution is shown. For illustration the convolution of a rectangular signal $x[k] = \\text{rect}_L[k]$ of length $L$ with a triangular signal $h[k] = \\Lambda_N[k]$ of length $N$ is considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2YXXV56P3v3UmgAygRg5YMxIDaCBZDelKF0noo0BOtVnO4kArF+lIlPg9arBokthX1qOEYqvY5tVqqIqC8FWOk6nGsIFZOLRgymKg5UUTeJihBGAE7QjK5nz/2GtgZZ/bsPbP3Xnvv+X6ua67Ze7381r1nzVrrnjX3+v0iM5EkSZJU8WtlByBJkiR1EhNkSZIkqYoJsiRJklTFBFmSJEmqYoIsSZIkVTFBliRJkqqYIEuacyLihIi4p4ntHR8RP4yIRyJi1STz74iIk4vX74yITzRr270kIjIinjXDdX8/IrY3OyZJc5MJsqRSFcnjaJFc/iQiPh0RB5QQw8mzaOK9wN9n5gGZubHWgpn5gcx8/Sy2JX41mc7Mb2bm0jJjktQ7TJAldYI/zswDgGOA5cDakuNp1DOA75UdRC0R0Vd2DJLULUyQJXWMzPwJMEglUQYgIvaNiAsj4q6I+GlEfDwi+ot5CyPiixExEhEPRMQ3I+LXinl73WEs7ky/b+I2I+IyYDHwL8Vd7HMniy0i3hARtxXbuTYiFhXTfwQcUbX+vrU+Y0S8OyI+U7xeUsT56uLz3R8Rf1W17K9FxHkR8aOI+FlEXB0RB1XN/+firvvPI+LfIuK5Ez7vxyLiyxHxC+APJonloIi4OCJ2RMSDEbGxat6kn7fqZ/vGoqzkwYj4aFTsW+yL36pa9uDiPwRPm67dCbHdEBGvr3r/moi4sXj9b8Xk7xQ/8z+ZWDYTEUcWbYxExPci4mUTfjYfjYgvRcTDEXFTRDyz1n6TNLeYIEvqGBFxKPBi4Laqyf8T+E0qSfOzgAHgXcW8twH3AAcDTwfeCWQj28zMVwF3UdzFzswPThLXicA64DTgEOBO4Mpi/WdOWP/RRrZf+D1gKXAS8K6IOLKY/hfAKuC/AouAB4GPVq33v4FnA08DNgOfndDuGcD7gScBN06y3cuA/YDnFm18eLrPW+WlwO8Ay4rlVhaffQNwetVypwHfyMz76mx3Wpn5wuLlsuJnflX1/IiYD/wL8NXic70Z+GxEVJdgnA68B3gKld+39zcah6TeZYIsqRNsjIiHgbuB+4DzASIigDcAf5mZD2Tmw8AHgFcW6+2ikmg9IzN3FXWoDSXIdfpT4FOZublIAtcCx0XEkia1/57MHM3M7wDfoZJ0AqwG/ioz7ym2+27g1IiYB5CZn8rMh6vmLYuIA6va/UJm/p/M3JOZv6zeYEQcQuWPkTdm5oPFz+8bDXzeCzJzJDPvAr7OE3f9L2fvBPmMYlq97TbDscABRYyPZeb1wBcnxLUhM2/OzN1U/rA4ZpJ2JM1RJsiSOsGqzHwScALwHGBhMf1gKnc4byn+VT4CfKWYDrCeyt2/r0bE7RFxXoviW0TlbicAmfkI8DMqd7Ob4SdVr/+TSnIHldrmz1d99m3AGPD0iOiLiAuK8ouHgDuKdRZWtXV3jW0eBjyQmQ9OMq+ezztVzNcD/RHxgoh4BpXE8/MNtNsMi4C7M3NP1bQ7qS9+STJBltQ5ijuYnwYuLCbdD4wCz83MBcXXgcUDfRR3T9+WmUcAfwy8NSJOKtb9TyrJ9bjfqLXpaULbQSVZBSAi9geeCgzX98lm7G7gxVWffUFm/npmDlO5M/ty4GTgQGDJeHhV69f6XHcDB0XEgknmzfjzFknp1VTu1p4BfLG4899ou7+g/v03WfyHjdejFxbXE78kgQmypM7zEeAPI+KYItn6J+DDVQ95DUTEyuL1SyPiWUUpxkNU7q6OFe3cCpxR3Gl9EZU63qn8lMqDdlO5HHhtRBxTPIT3AeCmzLxj5h+zLh8H3l/ciR1/4O3lxbwnAY9SuQO7XxFT3TLzXio1zP8QEU+JiPkRMV7bO9vPeznwJ1RKKi6fML3edm8FTomI/aLysOWfT5hfa5/dRCXBPrf4XCdQ+QOq4XpnSXOTCbKkjpKZO4FLgb8pJr2DShnFfxSlBF+j8kAbVB5Q+xrwCPAt4B8y84Zi3jlUkqIRKolarf6J1wF/XZQyvH2SmK4r4vkccC/wTJ6og26lvwOupVJC8jDwH8ALinmXUikbGAa+X8xr1Kuo1HH/Xyq132+B2X/ezBxPUBdRScLHpzfS7oeBx6gkwpfwqw8gvhu4pNhnp03Y/mPAy6jUWN8P/APwZ5n5f+v9DJLmtmjN8yySJElSd/IOsiRJklTFBFmSJEmqYoIsSZIkVTFBliRJkqrMKzuAWhYuXJhLliwpOwxJkiT1gFtuueX+zDx4uuU6OkFesmQJmzZtKjsMSZIk9YCIuHP6pSyxkCRJkvZigixJkiRVaWuJRUT8JfB6IIGtwGsz85ftjEGSesXGoWHWD25nx8goixb0s2blUlYtHyg7LEnqem1LkCNiAPgL4KjMHI2Iq6kMMfrpdsUgSWVqZkK7cWiYtRu2MrprDIDhkVHWbtgKMKs2Tbglqf0lFvOA/oiYB+wH7Gjz9iWpFOMJ7fDIKMkTCe3GoeEZtbd+cPvjyfG40V1jrB/c3hHxSVI3a1uCnJnDwIXAXcC9wM8z86sTl4uIsyJiU0Rs2rlzZ7vCk6SWanZCu2NktKHp02l2fJLUzdqWIEfEU4CXA4cDi4D9I+LMictl5kWZuSIzVxx88LTd1ElSV2h2QrtoQX9D06fT7PgkqZu1s8TiZODHmbkzM3cBG4DfbeP2JakhG4eGOf6C6zn8vC9x/AXXz6rcoNkJ7ZqVS+mf37fXtP75faxZuXRG7TU7vnHN/BlKUru0M0G+Czg2IvaLiABOAra1cfuSVLdm1+Q2O6FdtXyAdacczT59ldP4wIJ+1p1y9Iwfqmt2fGBds6Tu1c4a5JuAa4DNVLp4+zXgonZtX5Ia0eya3GYntONtLl+8gBccfhD/57wTZ91Ws+OzrllSt2prP8iZeT5wfju3KUkz0Yqa3FXLB7ji5rsAuGr1cTNup1WaHZ91zZK6lSPpSdIkWlWTO5f4M5TUrUyQJfWMZj4Q1oqa3LmmVXXNPvQnqdXaWmIhSa3S7JHlxtc595otPDa2hwFHlmtYs3+GrRg9UJImY4IsqSfUeiBspslTp9cMd4Nm/gxbsY8laTKWWEjqCT4Q1vvcx5LaxQRZUk/wgbDe5z6W1C4myJJ6gg/V9T73saR2sQZZUk/wobre5z6W1C4myJJKs3FomPWD29kxMsqiJiQ7PlTX+5q9j5v9OyipN5ggSyqFXXapbP4OSpqKNciSSlGryy6pHfwdlDQVE2RJpbDLLpXN30FJUzFBllQKu+xS2fwdlDQVE2RJpbDLLpXN30FJU/EhPUmlsMsulc3fQUlTMUGWVBq7ZVPZ/B2UNBlLLCRJkqQqJsiSJElSFUssJNXNUcek2jxGpN4wbYIcEQfV0c6ezBxpQjySOpSjjkm1eYxIvaOeO8g7iq+osUwfsLgpEUnqSLVGHfPiL3mMSL2kngR5W2Yur7VARAw1KR5JHcpRx6TaPEak3lHPQ3r19Htj3zhSj3PUMak2jxGpd0ybIGfmLwEi4n0T50VEX/UyknqXo45JtXmMSL2jkW7eBiLi9PE3EfE04GvND0lSJ1q1fIB1pxzNPn2V08bAgn7WnXK0tZVSwWNE6h2NdPO2GhiMiB8BCVwMvKORjUXEAuATwG8VbbwuM7/VSBuSyuOoY1JtHiNSb6inm7dLgc3AEHA2cDmwG1iVmbc1uL2/A76SmadGxD7Afg2uL0mSJLVUPSUWlxTLvY5KcrwEeBA4MyJOrXdDEfFk4IXAJwEy8zH7TpYkSVKnmfYOcmZeB1w3/j4i5gFHAcuAY4Fr6tzWEcBO4OKIWAbcApyTmb9oNGhJkiSpVRp5SA+AzNydmVsy87LMfHsDq84Dfhv4WNGv8i+A8yYuFBFnRcSmiNi0c+fORsOTJEmSZqWeGuTNmfnbs10GuAe4JzNvKt5fwyQJcmZeBFwEsGLFipwuPklT2zg0zPrB7ewYGWXRgn7WrFzqE/VSF/EYlspRTy8WR0bElhrzAzhwukYy8ycRcXdELM3M7cBJwPfrjFNSgzYODbN2w9bHh74dHhll7YatAF5gpS7gMSyVp54E+Tl1LDM2/SIAvBn4bNGDxe3Aa+tcT1KD1g9uf/zCOm501xjrB7d7cZW6gMewVJ56HtK7s1kby8xbgRXNak/S1HaMjDY0XVJn8RiWytPIQCEARMTVxcs7qfSPvLkomZDUQRYt6Gd4kgvpogX9JUQjqVEew1J5ZtKLxWmZeRrwceD3gG83PSpJs7Zm5VL65/ftNa1/fh9rVi4tKSJJjfAYlsozkzvIJwEvpTIK3jeBdzY7KEmzN16jeO41W3hsbA8DPgEvdRWPYak8DSfIwMXAV4AbgE2Z+fOmRiSpaVYtH+CKm+8C4KrVx5UcjaRGeQxL5Wg4Qc7MxRFxKPBfqAw3/ezMPL35oUmSJEntV3eCHBEnAn8KjADfBbYAX8nMR1sUmyRJktR2jdxB/gxwdrHO84BVwHOBZ7UgLkmSJKkUjSTIt2Xm54vX/9yKYCRJkqSyTdvNW0RcGhFvAb4VEW9rQ0ySJElSaerpB/mSYrnfAF4VEXdGxLUR8T8i4hWtDU+SJElqr3qGmr4OuG78fUTMA44ClgHPx3ILqWk2Dg2zfnA7O0ZGWWSfp5KazHOMVJ+ZdPO2m0oPFluaH440d20cGmbthq2M7hoDYHhklLUbtgJ4AZM0a55jpPo1PNS0pNZYP7j98QvXuNFdY6wf3F5SRJJ6iecYqX51JcgRcUbx/ZWtDUeau3aMjDY0XZIa4TlGql+9d5AHIuI04NBWBiPNZYsW9Dc0XZIa4TlGql893bydDxwEXA4cFBHvanlU0hy0ZuVS+uf37TWtf34fa1YuLSkiSb3Ec4xUv2kT5Mx8D/AAcCbwQGa+t+VRSXPQquUDrDvlaPbpqxyWAwv6WXfK0T48I6kpPMdI9au3F4sdmXllRJze0mikOW7V8gGuuPkuAK5afVzJ0UjqNZ5jpPrUVYOcmZ8tvl/R2nAkSZKkctnNmyRJklRl2hKLiDiojnb2ZOZIE+KRJEmSSlVPDfKO4itqLNMHLG5KRJIkSVKJ6kmQt2Xm8loLRMRQk+KRJEmSSlVPDfJxABHxvokzIqKvehlJkiSp29XTD/Ivi5cD40NOA0TE04CvTVhGkiRJ6mr19oMMsBoYjIjbgAQuBt7RkqgkSZKkktTTi8WlwGZgCDibypDTu4FVmXlboxssyjI2AcOZ+dJG15c6xcahYdYPbmfHyCiLFvSzZuVSR6SSNKd4HlSvqqcG+ZJiuddRSY6XAA8CZ0bEqTPY5jnAthmsJ3WMjUPDrN2wleGRURIYHhll7YatbBwaLjs0SWoLz4PqZfXUIF+XmR/KzFdn5jHAQuBtwI+AYxvZWEQcCrwE+MRMgpU6xfrB7YzuGttr2uiuMdYPbi8pIklqL8+D6mWN1CADkJm7gS3F12UNrv4R4FzgSVMtEBFnAWcBLF5s18rqTDtGRhuaLkm9xvOgetm0d5AjYnOTlnkpcF9m3lJrucy8KDNXZOaKgw8+eLpmpVIsWtDf0HRJ6jWeB9XL6qlBPjIittT42kql7GI6xwMvi4g7gCuBEyPiM7OIXSrNmpVL6Z/ft9e0/vl9rFm5tKSIJKm9PA+ql9VTYvGcOpYZm26BzFwLrAWIiBOAt2fmmXW0LXWc8ae0z71mC4+N7WHAp7clzTGeB9XLpk2QM/POdgQidZtVywe44ua7ALhqtYNJSpp7PA+qV9VTYlFTRPx+o+tk5g32gSxJkqRONOsEGXhFE9qQJEmSOkLD3bxFxLXAj6mMrnfLTNqQJEmSOlXdd5Aj4iMREZn5MuBDwEPAnwDPaFVwkiRJUrs1UmLxCHBtROxXPLj3C+DEzHxJa0KTJEmS2q/u8ojM/OuIOAP4RkQ8SiVBPq9lkUmSJEklqDtBjoiTgDdQSYwPAf48Mx1wXZIkST2lkRKLvwL+JjNPAE4FroqIE1sSlSRJklSSRkosTqx6vTUiXgx8DvjdVgQmNdvGoWHWD25nx8goixzxSZI6judpdYoZd9GWmfcWZRdSx9s4NMzaDVsZ3VUZFX14ZJS1G7YCePKVpA7geVqdZFYDhWTmaLMCkVpp/eD2x0+640Z3jbF+0DJ6SeoEnqfVSZoxkp7U8XaMTP633FTTJUnt5XlanaShBHn8oTwfzlO3WbSgv6HpkqT28jytTtLoHeQLJ3yXusKalUvpn9+317T++X2sWbm0pIgkSdU8T6uTzPQhvWhqFFKLjT/gce41W3hsbA8DPh0tSR3F87Q6yYx7sZC6zarlA1xx810AXLX6uJKjkSRN5HlancKH9CRJkqQqJsiSJElSlUYT5EeK7w83OxBJkiSpEzSUIGfmC6u/S5IkSb3GEgtJkiSpigmyJEmSVGXabt4i4qA62tmTmSNNiEd63MahYdYPbmfHyCiL7A9TktQgryOaqXr6Qd5RfNUaHKQPWNyUiCQqJ7W1G7YyumsMgOGRUdZu2ArgyU2SNC2vI5qNekostmXmEZl5+FRfwM9aHajmlvWD2x8/qY0b3TXG+sHtJUUkSeomXkc0G/UkyMcBRMT7Js6IiL7qZaRm2TEy2tB0SZKqeR3RbEybIGfmL4uXAxFxxvj0iHga8LUJy0wpIg6LiK9HxLaI+F5EnDPToNX7Fi3ob2i6JEnVvI5oNhrpxWI18IaIeH5E/A5wPXBhA+vvBt6WmUcCxwJnR8RRDayvOWTNyqX0z+/ba1r//D7WrFxaUkSSpG7idUSzUU8vFpcCm4Eh4GzgcirJ7qrMvK3eDWXmvcC9xeuHI2IbMAB8fwZxq8eNP0Bx7jVbeGxsDwM+fSxJaoDXEc1GPb1YXAIsA15XfF8CfBs4MyK+m5nXNLrRiFgCLAdummTeWcBZAIsX2zHGXLZq+QBX3HwXAFettsxdktQYryOaqWkT5My8Drhu/H1EzAOOopIsHws0lCBHxAHA54C3ZOZDk2zvIuAigBUrVmQjbUuSJEmzVc8d5L1k5m5gS/F1WSPrRsR8KsnxZzNzQ6PbliRJklpt2of0ImJzk5YJ4JNU+lX+UH3hSZIkSe1Vzx3kIyNiS435ARxYRzvHA68CtkbErcW0d2bml+tYV5IkSWqLehLk59SxzNh0C2TmjdQerlpdzjHvJUm9zmvd3FDPQ3p3AkTE16j0Y/ydlkelruOY95KkXue1bu5oZKCQc4EPR8TFEXFIqwJSd3LMe0lSr/NaN3fUnSBn5ubMPBH4IvCViDg/IhyvUYBj3kuSep/XurmjkTvI4z1RbAc+BrwZ+GFEvKoVgam7OOa9JKnXea2bO+pOkCPiRmAY+DCVIaJfA5wAPD8iLmpFcOoejnkvSep1XuvmjkYGCnkj8L3MnDi63ZsjYlsTY1IXcsx7SVKv81o3d9SdIGfmd2vMfkkTYlGXc8x7SVKv81o3NzRUgzyVzLy9Ge1IkiRJZWukxEI9xs7OJUkql9fizmSCPEfZ2bkkSeXyWty5mlJioe5jZ+eSJJXLa3HnMkGeo+zsXJKkcnkt7lwmyHOUnZ1LklQur8WdywR5jrKzc0mSyuW1uHP5kN4cZWfnkiSVy2tx5zJBnsPs7FySpHJ5Le5MJshdxL4SJUlSLeYKzWGC3CXsK1GSJNVirtA8PqTXJewrUZIk1WKu0DwmyF3CvhIlSVIt5grNY4LcJewrUZIk1WKu0DwmyC20cWiY4y+4nsPP+xLHX3A9G4eGZ9yWfSVKkqRaWpErNDOX6SY+pNcizS6Ut69ESZJUS7Nzhbn80J8JcovUKpSf6S+VfSVKkqRampkrtCKX6RYmyFWa2XeghfKSJKmbtSKX6ZZ+mttagxwRL4qI7RFxW0Sc185tT2f83wjDI6MkT/wbYaa1NhbKS5KkbtbsXKbZuVYrtS1Bjog+4KPAi4GjgNMj4qh2bX86ze470IfqJElSN2t2LtNN/TRHZrZnQxHHAe/OzJXF+7UAmbluqnVWrFiRmzZtakt8h5/3JRJYveULHPHzvf+SOfaIp86ozfsfeZQf7fwFmcm+8/o47KB+Fh6w76zi/P69DwFw1CFPnlU7tteZ7bWiTduzvTLba0Wbtmd7ZbbXijY7ub1m5jL/cfvPHn99+4ED/OPzXg5AAD++4CWzjrUeEXFLZq6Ybrl21iAPAHdXvb8HeMHEhSLiLOAsgMWLF7cnMir/LhiepKZm33l9kyxdn4UH7Msjj+4GYMlT959xO9X222fm8dhe57fXijZtz/bKbK8Vbdqe7ZXZXiva7OT2mpnL7Duvj0d3j/3K9E4sP23nHeRXACsz8/XF+1cBz8/MN0+1TjvvIE/sygQq/0ZYd8rRHVk8LkmS1E06IdfqxDvI9wCHVb0/FNjRxu3XNL5juuHJSkmSpG7TTblWO+8gzwN+AJwEDAPfBs7IzO9NtU477yBLkiSpt3XcHeTM3B0RbwIGgT7gU7WSY0mSJKkMbR0oJDO/DHy5nduUJEmSGtHWgUIkSZKkTte2GuSZiIidwJ0lbHohcH8J29Xk3B+dx33SWdwfncd90lncH52nrH3yjMw8eLqFOjpBLktEbKqngFvt4f7oPO6TzuL+6Dzuk87i/ug8nb5PLLGQJEmSqpggS5IkSVVMkCd3UdkBaC/uj87jPuks7o/O4z7pLO6PztPR+8QaZEmSJKmKd5AlSZKkKibIkiRJUhUT5CoR8aKI2B4Rt0XEeWXHI4iIOyJia0TcGhGbyo5nLoqIT0XEfRHx3appB0XEv0bED4vvTykzxrlkiv3x7ogYLo6TWyPij8qMcS6JiMMi4usRsS0ivhcR5xTTPUZKUmOfeJyUICJ+PSJujojvFPvjPcX0wyPipuIYuSoi9ik71mrWIBciog/4AfCHwD3At4HTM/P7pQY2x0XEHcCKzLSD95JExAuBR4BLM/O3imkfBB7IzAuKPyafkpnvKDPOuWKK/fFu4JHMvLDM2OaiiDgEOCQzN0fEk4BbgFXAa/AYKUWNfXIaHidtFxEB7J+Zj0TEfOBG4BzgrcCGzLwyIj4OfCczP1ZmrNW8g/yE5wO3ZebtmfkYcCXw8pJjkkqXmf8GPDBh8suBS4rXl1C5+KgNptgfKklm3puZm4vXDwPbgAE8RkpTY5+oBFnxSPF2fvGVwInANcX0jjtGTJCfMADcXfX+HjygOkECX42IWyLirLKD0eOenpn3QuViBDyt5HgEb4qILUUJhv/OL0FELAGWAzfhMdIRJuwT8DgpRUT0RcStwH3AvwI/AkYyc3exSMflXCbIT4hJpll/Ur7jM/O3gRcDZxf/Xpa0t48BzwSOAe4F/rbccOaeiDgA+Bzwlsx8qOx4NOk+8TgpSWaOZeYxwKFU/mN/5GSLtTeq2kyQn3APcFjV+0OBHSXFokJm7ii+3wd8nsqBpfL9tKjzG6/3u6/keOa0zPxpcQHaA/wTHidtVdRVfg74bGZuKCZ7jJRosn3icVK+zBwBbgCOBRZExLxiVsflXCbIT/g28Oziqcp9gFcC15Yc05wWEfsXD1gQEfsD/w34bu211CbXAq8uXr8a+EKJscx544lY4b/jcdI2xQNInwS2ZeaHqmZ5jJRkqn3icVKOiDg4IhYUr/uBk6nUhX8dOLVYrOOOEXuxqFJ0+fIRoA/4VGa+v+SQ5rSIOILKXWOAecDl7pP2i4grgBOAhcBPgfOBjcDVwGLgLuAVmemDY20wxf44gcq/jRO4A1g9Xv+q1oqI3wO+CWwF9hST30ml5tVjpAQ19snpeJy0XUQ8j8pDeH1UbsxenZnvLa7xVwIHAUPAmZn5aHmR7s0EWZIkSapiiYUkSZJUxQRZkiRJqmKCLEmSJFUxQZYkSZKqmCBLkiRJVUyQJUmSpComyJIkSVIVE2RJ6lIRcXJEXFZ2HJLUa0yQJal7LaMyApUkqYlMkCWpey0DhiJi34j4dER8ICKi7KAkqdvNKzsASdKMLQPuAwaBT2TmZ0qOR5J6QmRm2TFIkhoUEfOB+4E7gdWZ+a2SQ5KknmGJhSR1p6OAbwO7gbGSY5GknmKCLEndaRnw78ArgYsj4uklxyNJPcMEWZK60zLgu5n5A+AdwNVF2YUkaZasQZYkSZKqeAdZkiRJqmKCLEmSJFUxQZYkSZKqmCBLkiRJVUyQJUmSpComyJIkSVIVE2RJkiSpigmyJEmSVMUEWZIkSapigixJkiRVMUGWJEmSqpggS5IkSVVMkCWpEBEnRMQ9TWzv+Ij4YUQ8EhGrJpm/NCKGIuLhiPiLZm2300XEayLixlms//GI+JtmxiRJ1UyQJXWkiLgjIkaL5PInEfHpiDighBhOnkUT7wX+PjMPyMyNk8w/F7ghM5+Umf/fTDcSETdExOtnHGUHmyyZzsw3Zub/KCsmSb3PBFlSJ/vjzDwAOAZYDqwtOZ5GPQP43izmS5JKYIIsqeNl5k+AQSqJMgARsW9EXBgRd0XET4t/u/cX8xZGxBcjYiQiHoiIb0bErxXzMiKeVdXOpyPifRO3GRGXAYuBfynuYp87WWwR8YaIuK3YzrURsaiY/iPgiKr1952w3vXAHwB/X8z/zYh4SVFy8VBE3B0R765a/tcj4jMR8bPic307Ip4eEe8Hfr+qnb+fIs7fi4h/L9a9OyJeU0w/MCIujYidEXFnRPx11c/qNRFxY/FzfjAifhwRLy7mvTIiNk3Yxl9GxLXTtTthnSXFPplXNe2GiHh9RBwJfBw4rvhsI5Pts6n2QTEvI+KNRanLgxHx0YiIyX5GkjTOBFlSx4uIQ4EXA7dVTf6fwG9SSZqfBQwA7yrmvQ24BzgYeDrwTiAb2WZmvgq4i+IudmZ+cJK4TgTWAacBhwB3AlcW6z9zwvqPTmj/ROCbwJuK+T8AfgH8GbAAeAnw/1TVLr8aOBA4DHgq8EZgNDP/akI7b5okzsXA/wb+V/EzOQa4tZj9v4pRnCy/AAAayElEQVR2jwD+a7H911at/gJgO7AQ+CDwySLBvBZYGhHPrlr2DODyOtudVmZuKz7nt4rPtmCSzzblPqjyUuB3gGXFcisbiUPS3GOCLKmTbYyIh4G7gfuA8wGKBO0NwF9m5gOZ+TDwAeCVxXq7qCRLz8jMXZn5zcxsKEGu058Cn8rMzUUCvJbK3c4lM2ksM2/IzK2ZuScztwBXUEkuofKZngo8KzPHMvOWzHyogTi/lplXFD+Pn2XmrRHRB/wJsDYzH87MO4C/BV5Vte6dmflPmTkGXELl5/r0zPxP4AvA6QBFovwc4No6222WevbBBZk5kpl3AV+n6j8RkjQZE2RJnWxVZj4JOIFK8rWwmH4wsB9wS1EyMAJ8pZgOsJ7K3eavRsTtEXFei+JbROWOJQCZ+QjwMyp3sxsWES+IiK8XZQk/p3L3dPwzX0alzOTKiNgRER+MiPl1Nn0Y8KNJpi8E9qn+DMXr6vh/Mv6iSIoBxh+WvJwiQaZy93hjsUw97TZLPfvgJ1Wv/5Mn4pekSZkgS+p4mfkN4NPAhcWk+4FR4LmZuaD4OrB4oI/iruXbMvMI4I+Bt0bEScW6/0kluR73G7U2PU1oO6g8aAdAROxP5S7vcH2f7FdcTqV04bDMPJBK/W0AFHd+35OZRwG/S6Vs4M/qjPNu4JmTTL+fyp3pZ1RNW9xA/F8FFkbEMVQS5fHyikba/UXxfap90u59IEkmyJK6xkeAP4yIYzJzD/BPwIcj4mkAETEQESuL1y+NiGcVpRgPAWPFF1Rqb8+IiL6IeBFPlDBM5qdUamincjnw2og4pngI7wPATUVJwUw8CXggM38ZEc+ncleW4jP9QUQcXZQvPEQlAR3/TNPF+Vng5Ig4LSLmRcRTi5/jGHA18P6IeFJEPAN4K/CZeoLNzN3ANVTu2B8E/Gsxve52M3MnlWT2zGKfvI69k/mfAodGxD5ThNHsfSBJJsiSukORSF0KjA8Q8Q4qZRT/EREPAV8Dlhbznl28fwT4FvAPmXlDMe8cKneVR6jUr07WP/G4dcBfF2Ucb58kpuuKeD4H3EslsXvlxOUa8P8C7y3qrt9FJckc9xtUktGHgG3AN3gi4fw74NSil4Zf6U+5qL39IyoPLz5A5Y+EZcXsN1O5i3s7cCOVhPNTDcR8OXAy8M9FwjyukXbfAKyhUhrxXODfq+ZdT6UrvJ9ExP2TfLZm7wNJIlrz3IokSZLUnbyDLEmSJFUxQZYkSZKqmCBLkiRJVUyQJUmSpCrzyg6gloULF+aSJUvKDkOSJEk94JZbbrk/Mw+ebrmOTpCXLFnCpk2byg5DkiRJPSAi7px+KUssJEmSpL2YIEuSJElV2lpiERF/CbweSGAr8NrM/GU7Y5CkXrFxaJj1g9vZMTLKogX9rFm5lFXLB8oOS5K6XtsS5IgYAP4COCozRyPiairDgX66XTFIUpmamdBuHBpm7YatjO4aA2B4ZJS1G7YCzKpNE25Jan+JxTygPyLmAfsBO9q8fUkqxXhCOzwySvJEQrtxaHhG7a0f3P54cjxudNcY6we3d0R8ktTN2pYgZ+YwcCFwF3Av8PPM/OrE5SLirIjYFBGbdu7c2a7wJKmlmp3Q7hgZbWj6dJodnyR1s3aWWDwFeDlwODAC/HNEnJmZn6leLjMvAi4CWLFiRbYrPkmaqJklB81OaBct6Gd4knUXLeifUXvNjm+cZRuSulE7SyxOBn6cmTszcxewAfjdNm5fkurW7JKDqRLXmSa0a1YupX9+317T+uf3sWbl0hm11+z4wLINSd2rnQnyXcCxEbFfRARwErCtjduXpLo1u+Sg2QntquUDrDvlaPbpq5zGBxb0s+6Uo2d8d7bZ8YFlG5K6V9tKLDLzpoi4BtgM7AaGKEopJKnTNLvkYDxxPfeaLTw2toeBJpQbrFo+wBU33wXAVauPm3E7rYqvVWUbktRqbe0HOTPPB85v5zYlaSaaXeMLzU1oW6HZ8bXiZyhJ7eBIepJ6xsahYY6/4HoOP+9LHH/B9bOqdW1FycFc04qfYTP3sSRNpa13kCWpVZo9cEYrSg7mmmb/DFsxOIokTcYEWVJPqPVA2EyTp04viegGzfwZtmIfS9JkLLGQ1BN8IKz3uY8ltYsJsqSe0Ip+fNVZ3MeS2sUEWVJP8KG63uc+ltQu1iBL6gk+VNf73MeS2sUEWVLP8KG63uc+ltQOJsiSSrNxaJj1g9vZMTLKIu8GqgT+DkqajAmypFLYp63K5u+gpKn4kJ6kUtTq01ZqB38HJU3FBFlSKezTVmXzd1DSVEyQJZXCPm1VNn8HJU3FBFlSKezTVmXzd1DSVHxIT1Ip7NNWZfN3UNJUTJAllcY+bVU2fwclTcYSC0mSJKmKCbIkSZJUxRILSXVz1DGpNo8RqTdMmyBHxEF1tLMnM0eaEI+kDuWoY1JtHiNS76jnDvKO4itqLNMHLG5KRJI6Uq1Rx7z4Sx4jUi+pJ0HelpnLay0QEUNNikdSh3LUMak2jxGpd9TzkF49/d7YN47U4xx1TKrNY0TqHdMmyJn5S4CIeN/EeRHRV72MpN7lqGNSbR4jUu9opJu3gYg4ffxNRDwN+FrzQ5LUiVYtH2DdKUezT1/ltDGwoJ91pxxtbaVU8BiRekcj3bytBgYj4kdAAhcD72hkYxGxAPgE8FtFG6/LzG810oak8jjqmFSbx4jUG+rp5u1SYDMwBJwNXA7sBlZl5m0Nbu/vgK9k5qkRsQ+wX4PrS5IkSS1VT4nFJcVyr6OSHC8BHgTOjIhT691QRDwZeCHwSYDMfMy+kyVJktRppr2DnJnXAdeNv4+IecBRwDLgWOCaOrd1BLATuDgilgG3AOdk5i+qF4qIs4CzABYvtmtlSZIktVfDQ01n5m5gS/F1WYPb+m3gzZl5U0T8HXAe8DcT2r8IuAhgxYoV2Wh8kp7gsLdSd/MYlsoxbYlFRGxuxjLAPcA9mXlT8f4aKgmzpBYYH/Z2eGSU5IlhbzcODZcdmqQ6eAxL5amnBvnIiNhS42srsHC6RjLzJ8DdETHeIeRJwPdnEbukGmoNeyup83kMS+Wpp8TiOXUsMzb9IgC8Gfhs0YPF7cBr61xPUoMc9lbqbh7DUnnqeUjvzmZtLDNvBVY0qz1JU1u0oJ/hSS6kDnsrdQePYak8jYykB0BEXF18rY+I06tKJiR1EIe9lbqbx7BUnpn0YnEaQEQ8E3gr8I/Ak5scl6RZGn/S/dxrtvDY2B4GfAJe6ioew1J5Gk6QI+Ik4KVURsH7JvDOZgclqTkc9lbqbh7DUjkaLrEALgb2B74B3JyZP29uSJIkSVJ5ZlJisTgiDgX+C5Xhpp+dmac3PzRJkiSp/epOkCPiROBPgRHgu1RG0vtKZj7aotgkSZKktmvkDvJngLOLdZ4HrAKeCzyrBXFJkiRJpWgkQb4tMz9fvP7nVgQjSZIklW3ah/Qi4tKIeAvwrYh4WxtikiRJkkpTzx3kS4BlwG8AKyPiL4DvFF9bMtO7yVKTbBwaZv3gdnaMjLLIPk8lNZnnGKk+9Qw1fR1w3fj7iJgHHEUlaX4+lltITbFxaJi1G7YyumsMgOGRUdZu2ArgBUzSrHmOkerXcD/Imbk7M7dk5mWZuaYVQUlz0frB7Y9fuMaN7hpj/eD2kiKS1Es8x0j1m8lAIZJaYMfIaEPTJakRnmOk+tWVIEfEGcX3V7Y2HGnuWrSgv6HpktQIzzFS/eq9gzwQEacBh7YyGGkuW7NyKf3z+/aa1j+/jzUrl5YUkaRe4jlGql893bydDxwEXA4cFBHvanlU0hy0avkA6045mn36KoflwIJ+1p1ytA/PSGoKzzFS/erpxeI9EbEGOBM4NDMvbH1Y0ty0avkAV9x8FwBXrT6u5Ggk9RrPMVJ96i2x2JGZVwLDrQxGkiRJKltdCXJmfrb4fkVrw5EkSZLKZTdvkiRJUpVpa5Aj4qA62tmTmSNNiEeSJEkq1bQJMrCj+Ioay/QBi5sSkSRJklSiehLkbZm5vNYCETHUpHgkSZKkUtVTg3wcQES8b+KMiOirXkaSJEnqdtMmyJn5y+LlwPiQ0wAR8TTgaxOWkSRJkrpaPSUW41YDgxFxG5DAxcA7Gt1gcdd5EzCcmS9tdH2pU2wcGmb94HZ2jIyyaEE/a1YudUQqSXOK50H1qnp6sbgU2AwMAWdTGXJ6N7AqM2+bwTbPAbYBT57BulJH2Dg0zNoNWxndNQbA8MgoazdsBfDiIGlO8DyoXlZPDfIlxXKvo5IcLwEeBM6MiFMb2VhEHAq8BPhEY2FKnWX94PbHLwrjRneNsX5we0kRSVJ7eR5UL5v2DnJmXgdcN/4+IuYBRwHLgGOBaxrY3keAc4EnTbVARJwFnAWweLE9x6kz7RgZbWi6JPUaz4PqZQ2PpJeZuzNzS2Zelplvr3e9iHgpcF9m3jJN+xdl5orMXHHwwQc3Gp7UFosW9Dc0XZJ6jedB9bJpE+SI2NyMZYDjgZdFxB3AlcCJEfGZOtaTOs6alUvpn9+317T++X2sWbm0pIgkqb08D6qX1dOLxZERsaXG/AAOnK6RzFwLrAWIiBOAt2fmmfUEKXWa8QdQzr1mC4+N7WHAp7clzTGeB9XL6kmQn1PHMmPTLyL1llXLB7ji5rsAuGq1Y+VImns8D6pX1fOQ3p3N3mhm3gDc0Ox2JUmSpNlq+CG9iSLi95sRiCRJktQJZp0gA69oQhuSJElSR2hkqGkAIuJa4MdURte7ZSZtSJIkSZ2q7jvIEfGRiIjMfBnwIeAh4E+AZ7QqOEmSJKndGimxeAS4NiL2Kx7c+wVwYma+pDWhSZIkSe1Xd3lEZv51RJwBfCMiHqWSIJ/XssgkSZKkEtSdIEfEScAbqCTGhwB/npnbWxWYJEmSVIZGHrD7K+BvMvPGiDgauCoi3pqZ17coNqmpNg4Ns35wOztGRlnkiE+S1HE8T6tTNFJicWLV660R8WLgc8DvtiIwqZk2Dg2zdsNWRndVBn0cHhll7YatAJ58JakDeJ5WJ5lxP8iZeS9wUhNjkVpm/eD2x0+640Z3jbF+0CohSeoEnqfVSWY1UEhmjjYrEKmVdoxM/qs61XRJUnt5nlYnacZIelLHW7Sgv6HpkqT28jytTtJQghwRJ1Z/l7rFmpVL6Z/ft9e0/vl9rFm5tKSIJEnVPE+rkzR6B/nCCd+lrrBq+QDrTjmaffoqv/IDC/pZd8rRPvghSR3C87Q6SSPdvFWLpkYhtcGq5QNccfNdAFy1+riSo5EkTeR5Wp3CGmRJkiSpigmyJEmSVMUEWZIkSarSaIL8SPH94WYHIkmSJHWChhLkzHxh9XdJkiSp11hiIUmSJFWZaTdvUsttHBpm/eB2doyMsmhBP2tWLrU/TElS3byOaKamTZAj4qA62tmTmSNNiEcCKie1tRu2MrprDIDhkVHWbtgK4MlNkjQtryOajXruIO8ovmoNDtIHLG5KRBKwfnD74ye1caO7xlg/uN0TmyRpWl5HNBv1JMjbMnN5rQUiYqhJ8UgA7BgZbWi6JEnVvI5oNup5SO84gIh438QZEdFXvUwtEXFYRHw9IrZFxPci4pzGQtVcsmhBf0PTJUmq5nVEszFtgpyZvyxeDkTEGePTI+JpwNcmLFPLbuBtmXkkcCxwdkQc1XjImgvWrFxK//y+vab1z+9jzcqlJUUkSeomXkc0G430YrEaGIyI24AELgbeUe/KmXkvcG/x+uGI2AYMAN9vIAbNEeP1Yedes4XHxvYw4NPHkqQGeB3RbNTTi8WlwGZgCDgbuJzK3eBVmXnbTDYaEUuA5cBNM1lfc8Oq5QNccfNdAFy1etoqHkmS9uJ1RDNVTw3yJcVyr6OSHC8BHgTOjIhTG91gRBwAfA54S2Y+NMn8syJiU0Rs2rlzZ6PNS5IkSbMy7R3kzLwOuG78fUTMA44CllGpJb6m3o1FxHwqyfFnM3PDFNu7CLgIYMWKFVlv25IkSVIzNDySXmbuBrYUX5fVu15EBPBJKt3GfajR7UqSJEntMG2JRURsbsYywPHAq4ATI+LW4uuP6lhPkiRJapt67iAfGRFbaswP4MDpGsnMG6k9Gp+6nGPeS5J6nde6uaGeBPk5dSwzNv0i6mWOeS9J6nVe6+aOegYKuTMz76RSP7xg/P2Er3taH6o6Wa0x7yVJ6gVe6+aOerp5G3cu8OGIuDgiDmlVQOpOjnkvSep1XuvmjroT5MzcnJknAl8EvhIR50eEA5oLcMx7SVLv81o3dzRyB3m8q7btwMeANwM/jIhXtSIwdRfHvJck9TqvdXNH3QlyRNwIDAMfBgaA1wAnAM+PiItaEZy6x6rlA6w75Wj26av8Sg0s6GfdKUf70IIkqWd4rZs7Ghko5I3A9zJz4uh2b46IbU2MSV3KMe8lSb3Oa93cUHeCnJnfrTH7JU2IRZIkSSpdQzXIU8nM25vRjiRJklS2Rkos1GMcDUiSpHJ5Le5MJshzlKMBSZJULq/FnaspJRbqPo4GJElSubwWdy4T5DnK0YAkSSqX1+LOZYI8RzkakCRJ5fJa3LlMkOcoRwOSJKlcXos7lw/pzVHjxf/nXrOFx8b2MOCTs5IktZXX4s5lgtxFmt0VjKMBSZJUrmZfi+02rjlMkLuEXcFIkqRazBWaxxrkLmFXMJIkqRZzheYxQe4SdgUjSZJqMVdoHhPkLmFXMJIkqRZzheYxQe4SdgUjSZJqMVdoHh/Sa6FmPklqVzCSJKmWVuQKc7VXDBPkFmnFk6R2yyZJkmppZq4wl3vFsMSiRXySVJIkdbO5nMt4B7lKM/+N4JOkkiSpm7Uil+mWko223kGOiBdFxPaIuC0izmvntqcz/m+E4ZFRkif+jbBxaHhG7fkkqSRJ6mbNzmWanWu1UtsS5IjoAz4KvBg4Cjg9Io5q1/an0+x/I/gkqSRJ6mbNzmW6qWQjMrM9G4o4Dnh3Zq4s3q8FyMx1U62zYsWK3LRpU1viO/y8L5HA6i1f4Iif7/2XzLFHPHVGbd7/yKP8aOcvyEz2ndfHYQf1s/CAfWcV5/fvfQiAow558qzasb3ObK8Vbdqe7ZXZXivatD3bK7O9VrTZye01M5f5j9t/9vjr2w8c4B+f93IAAvjxBS+Zdaz1iIhbMnPFdMu1swZ5ALi76v09wAsmLhQRZwFnASxevLg9kVH5d8HwJDU1+87rm2Tp+iw8YF8eeXQ3AEueuv+M26m23z4zj8f2Or+9VrRpe7ZXZnutaNP2bK/M9lrRZie318xcZt95fTy6e+xXpndi+Wk77yC/AliZma8v3r8KeH5mvnmqddp5B3liVyZQ+TfCulOO7sjicUmSpG7SCblWJ95Bvgc4rOr9ocCONm6/pvEd0w1PVkqSJHWbbsq12nkHeR7wA+AkYBj4NnBGZn5vqnXaeQdZkiRJva3j7iBn5u6IeBMwCPQBn6qVHEuSJEllaOtAIZn5ZeDL7dymJEmS1AiHmpYkSZKqtK0GeSYiYidwZwmbXgjcX8J2NTn3R+dxn3QW90fncZ90FvdH5ylrnzwjMw+ebqGOTpDLEhGb6ingVnu4PzqP+6SzuD86j/uks7g/Ok+n7xNLLCRJkqQqJsiSJElSFRPkyV1UdgDai/uj87hPOov7o/O4TzqL+6PzdPQ+sQZZkiRJquIdZEmSJKmKCbIkSZJUxQS5SkS8KCK2R8RtEXFe2fEIIuKOiNgaEbdGxKay45mLIuJTEXFfRHy3atpBEfGvEfHD4vtTyoxxLplif7w7IoaL4+TWiPijMmOcSyLisIj4ekRsi4jvRcQ5xXSPkZLU2CceJyWIiF+PiJsj4jvF/nhPMf3wiLipOEauioh9yo61mjXIhYjoA34A/CFwD/Bt4PTM/H6pgc1xEXEHsCIz7eC9JBHxQuAR4NLM/K1i2geBBzLzguKPyadk5jvKjHOumGJ/vBt4JDMvLDO2uSgiDgEOyczNEfEk4BZgFfAaPEZKUWOfnIbHSdtFRAD7Z+YjETEfuBE4B3grsCEzr4yIjwPfycyPlRlrNe8gP+H5wG2ZeXtmPgZcCby85Jik0mXmvwEPTJj8cuCS4vUlVC4+aoMp9odKkpn3Zubm4vXDwDZgAI+R0tTYJypBVjxSvJ1ffCVwInBNMb3jjhET5CcMAHdXvb8HD6hOkMBXI+KWiDir7GD0uKdn5r1QuRgBTys5HsGbImJLUYLhv/NLEBFLgOXATXiMdIQJ+wQ8TkoREX0RcStwH/CvwI+AkczcXSzScTmXCfITYpJp1p+U7/jM/G3gxcDZxb+XJe3tY8AzgWOAe4G/LTecuSciDgA+B7wlMx8qOx5Nuk88TkqSmWOZeQxwKJX/2B852WLtjao2E+Qn3AMcVvX+UGBHSbGokJk7iu/3AZ+ncmCpfD8t6vzG6/3uKzmeOS0zf1pcgPYA/4THSVsVdZWfAz6bmRuKyR4jJZpsn3iclC8zR4AbgGOBBRExr5jVcTmXCfITvg08u3iqch/glcC1Jcc0p0XE/sUDFkTE/sB/A75bey21ybXAq4vXrwa+UGIsc954Ilb473ictE3xANIngW2Z+aGqWR4jJZlqn3iclCMiDo6IBcXrfuBkKnXhXwdOLRbruGPEXiyqFF2+fAToAz6Vme8vOaQ5LSKOoHLXGGAecLn7pP0i4grgBGAh8FPgfGAjcDWwGLgLeEVm+uBYG0yxP06g8m/jBO4AVo/Xv6q1IuL3gG8CW4E9xeR3Uql59RgpQY19cjoeJ20XEc+j8hBeH5Ubs1dn5nuLa/yVwEHAEHBmZj5aXqR7M0GWJEmSqlhiIUmSJFUxQZYkSZKqmCBLkiRJVUyQJUmSpComyJIkSVIVE2RJkiSpigmyJEmSVMUEWZK6VEScHBGXlR2HJPUaE2RJ6l7LqIxAJUlqIhNkSepey4ChiNg3Ij4dER+IiCg7KEnqdvPKDkCSNGPLgPuAQeATmfmZkuORpJ4QmVl2DJKkBkXEfOB+4E5gdWZ+q+SQJKlnWGIhSd3pKODbwG5grORYJKmnmCBLUndaBvw78Erg4oh4esnxSFLPMEGWpO60DPhuZv4AeAdwdVF2IUmaJWuQJUmSpCreQZYkSZKqmCBLkiRJVUyQJUmSpComyJIkSVIVE2RJkiSpigmyJEmSVMUEWZIkSary/wOyBgaem8K9IgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 16  # length of signal x[k]\n",
    "N = 16  # length of signal h[k]\n",
    "M = N+L-1\n",
    "\n",
    "# generate signals\n",
    "x = np.ones(L)\n",
    "h = sig.triang(N)\n",
    "\n",
    "# linear convolution\n",
    "y1 = np.convolve(x, h, 'full')\n",
    "# fast convolution\n",
    "y2 = np.fft.ifft(np.fft.fft(x, M) * np.fft.fft(h, M))\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.subplot(211)\n",
    "plt.stem(y1)\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$y[k] = x_L[k] * h_N[k]$')\n",
    "plt.title('Result of linear convolution')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.stem(y1)\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$y[k] = x_L[k] * h_N[k]$')\n",
    "plt.title('Result of fast convolution')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Numerical complexity\n",
    "\n",
    "It was already argued that the numerical complexity of the fast convolution is considerably lower due to the usage of the FFT. The gain with respect to the convolution is evaluated in the following. In order to measure the execution times for both algorithms the `timeit` module is used. The algorithms are evaluated for the convolution of two random signals $x_L[k]$ and $h_N[k]$ of length $L=N=2^n$ for $n=0, 1, \\dots, 16$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4sAAAJgCAYAAAAj5/aWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XuYXWV99//3xwRCCEoUECNRoj+x9RBEBesRp2jVQiu11bY2UrD1ibYea6zFPm09tD7FAxSt9ddGRalFURGPeKK2G+oBbMBAwOAJo5wUAQcJVnTi9/ljrXnczJ5kBid71hzer+uaK3vd67C/9569rj2f3PdaO1WFJEmSJEn97tR1AZIkSZKkucewKEmSJEkaYFiUJEmSJA0wLEqSJEmSBhgWJUmSJEkDDIuSJEmSpAGGRUnSTiVZl+Qzs/yc25I8cTafs0tJekme0z6e9dd7Qi1/l+SGJN/tqoYuzOQ9l+TeSbYnWbK765KkrhkWJWkWJPmDJJvaPyqvS/LJJI/tuq6pVNUZVfWkruuYrvkeNIf1eidZk6SSLJ1k3cYk65PcC9gAPLCq7jGD5xpJcvVM6p3LJr7Hquo7VbVPVe3osi5JGgbDoiQNWZKXAqcC/wc4ELg38Fbg2C7rmspkwULdGeLv4ynAJ4CDgRur6vohPY8kaZ4xLErSECXZF3gN8PyqOruqbq2qn1bVx6rqz9ttliU5Ncm17c+pSZa160aSXJ3k5Umub0clfyvJ0Um+luSmJH/Z93yvSnJWkvcluSXJxUke0rf+xCTfbNd9JcnT+tadkOTzSf4hyU3Aq9q2z7Xr0667PsnNSS5N8uDxfib51yTfT/LtJH+V5E59x/1ckjcm+UGSbyX59SleuiPa+n6Q5J1J9uqr8zeSbE4ymuQLSQ5t299NE8Q/1o7gvjzJ6Uk2tOsPakfX/rRdvl/7+mVXx23X3TPJB9v+fSvJiya85u9v+39LksuTHL6L98SvJbmifQ3fAmTC7+BzfcuV5PlJvg58vW375STntrV/Ncnv9m2/PMnJ7e/g5vZ1Xw6c324y2r42j2q3PxQYBX4ZOBe4Z7v+Xe36DyT5bnus85M8qO+5jm5/R7ckuSbJy5KsAD7Zd5ztSe45yWuwszpJ8tT2NRxNM0X3AX37bWuf59J2v/eNvzeSbE3yG33bLk0zpfZhUx13Qm3vSvJ3fcv/b6R0J++x243atu+Vj7a/n28k+V99x7pD7xVJ6pphUZKG61HAXsCHdrHN/wYeCRwGPAR4BPBXfevv0R7jIOBvgLcBzwIeDjwO+Jsk9+3b/ljgA8DdgPcAH06yR7vum+0++wKvBv4tyaq+fX8FuBK4O/DaCXU+CTgSuD+wEvg94MZ23T+2x7wv8HjgD4FnTzjuV4H9gdcD7xgPaTuxDngy8P+1z/dXAO0f/qcBzwX2A/4F+GiSZVV1HPAd4DfbaYGvB84DRtpjPr7t2+Pb5SOB/6qq2tVx04TejwGX0PwOngC8JMmT++p9KnBm+7p8FHjLZJ1Ksj/wwbY/+9P8Ph6zi9cB4LdoXr8HtmHsXJrf692BZwJv7Qtxb6R5Xzya5vf/cuBnbV8BVravzRfb5aOBc6rq34FfB65t15/Qrv8kcEj7XBcDZ/TV9Q7guVV1Z+DBwH9U1a0TjrNPVV07SZ8mrTPJ/YH3Ai8BDqAZ8fxYkj379v1dmtHQ+wCHAuO1vrd9PcY9Gbihqi6e5nGntJP32ETvBa4G7gk8Hfg/SZ7Qt35a7xVJmgsMi5I0XPvR/ME6tott1gGvqarrq+r7NCHuuL71PwVeW1U/pfkjc3/gTVV1S1VdDlxO80fzuIuq6qx2+1NoguYjAarqA1V1bVX9rKreRzNa9Yi+fa+tqn+sqrGq+p8Jdf4UuDPNKFSqamtVXZfmxh6/B7yirWkbcPKEPny7qt7WXtd1OrCKZkruzrylqq6qqptoQut4CPhfwL9U1YVVtaOqTgduG+/fJM4DHtcGviNpgup4OHt8u36q4x4BHFBVr6mqn1TVlTSB/ff7nudzVfWJtn/vpgn9kzka+Erf7+dUYKqbyfx9Vd3U/j5+A9hWVe9sf0cX04TPp7d9/CPgxVV1TduPL1TVbbs49jE0wWlSVXVa+zu9DXgV8JA0o+XQvB8emOQuVfWDtpYpTVHn79GE13Pb1+eNwHKaUDnuze17+CaaEH9Y2/4e4KlJ9m6X/6BtY5rHnbE0130+FviLqvpxVW0G3s7tz4XpvlckqXOGRUkarhuB/bPr683uCXy7b/nbbdv/O0bfzTPGA9z3+tb/D7BP3/JV4w+q6mf8fJSDJH+Yn0+1HKUZEdp/sn0nqqr/oBkF+Sfge2lujHKXdv89J+nDQX3L3+07zo/ah/01T9RfR//rcTCwYbz+tg/34vavV3/N3wS20wSKxwEfB65N8kvcPizu6rgH00yr7F/3l9w+7PYHvh8Be+3kd35Pbv/7KXbxmk/yWhwM/MqEWtbRjD7vT/MfA9+c4ngAJFlJE/y/sJP1S5KclGba8g+Bbe2q8ffL79CE328nOW98aus07KrO250L7fv3KnbyXqJ5rfdpt/0GsBX4zTYwPpWfh8XpHHd3uCdwU1Xd0te203OBXb9XJKlzhkVJGq4vAj+mmUq4M9fShIBx927bflH3Gn/QjuKspglIB9OMiL0A2K+qVgKX0XfNHFC7OnBVvbmqHg48iGZ66J8DN9CMMk3swzW7ow/c/vW4imaUdWXfz95V9d5d1H8ezXTAPavqmnb5D4G7ApuncdyrgG9NWHfnqjr6F+jXddz+95MJfZ1Mf5+uAs6bUMs+VfUnNL+HH9NM3d3VMcY9GfjsLu7i+Qc0U5qfSDPFeM142QBV9d9VdSzNFNUPA+/fxXP121WdtzsX+l6f6b6XxqeiHkszgvuNX+C4twJ79y1PvDPsrvp3LXC3JHfua5vpuSBJnTEsStIQVdXNNNcZ/lOaG9PsnWSPJL+eZPx6p/cCf5XkgPaatr8B/m0GT/vwJL/djla8hGY65QXACpo/dL8PkOTZNCOL05LkiCS/0l7/eCvNH/w72rDxfuC1Se7chtKXzrAPz0+yOsndaEbx3te2vw14XltHkqxIckzfH+ffo7lust95NAF5/CYvPeCFNNMBx4PSro77JeCHSf4izY1ZliR5cJIjfoF+nQM8qO/38yIGw8iufBy4f5Lj2vfRHu3v5QHtaNlpwCntTVaWJHlUmpslfZ/m2sX+12aXU1BpphzfRjM6vjfN3XwBSLJnmu+E3Led1vlDYPy1/B6wX9901duZos73A8ckeUL7PtvQ1jDp6OckzqS5tvZP+PmoInfwuJuBo5PcLck9aM6hfpO9x8b7dlV7zL9PsleaGwj9Mbe/1lOS5g3DoiQNWVWdQhOe/ormj/araMLLh9tN/g7YBFwKbKG5kcjfDR5p2j5Cc43WD2iulfrtau7A+hWaawm/SPMH71rg83fguHehCVU/oJladyPNtV/QhK9baW4g8zmaP9RPm0Ef3gN8pj3elbSvR1Vtorm+8C1tHd/g5zc4Afh7muA9muRlbdt5NMFnPCx+jib8jC/v8rhtoPxNmqms36IZGXs7zWjbHVJVNwDPAE6ief0O4Q78DtrpjU+iuV7yWpopja8DlrWbvIzmPfTfwE3tuju1U39fC3y+fW0eBfwa8KldPN2/0vyerwG+QvMfDv2OA7a1U1SfR3PTJarqCpr/ALmyfa7JpgjvrM6vtsf5R5rX+Tdpbibzk6lfHaiq62je34/m5//BwB087rtpbma0jeY9+L4J6yd7j/V7Js0o7LU0N7Z6ZVWdO536JWmuSXO5hCRpIUjyKuB+VfWsrmvR3JXkETQ3EXrElBtLkhYtRxYlSVqcXtl1AZKkuc27b0mStMhU1Ze6rkGSNPc5DVWSJEmSNMBpqJIkSZKkAYtuGur+++9fa9as6bqMGbn11ltZsWJF12VIc5rniTQ1zxNpap4nWoguuuiiG6rqgKm2W3Rhcc2aNWzatKnrMmak1+sxMjLSdRnSnOZ5Ik3N80SamueJFqIk357Odk5DlSRJkiQNMCxKkiRJkgYYFiVJkiRJAwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjTAsChJkiRJGmBYlCRJkiQNMCxKkiRJkgYYFiVJkiRJAwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjTAsChJkiRJGmBYlCRJkiQNMCxKkiRJkgYYFiVJkiRJAwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpwNKuC5Ck3e4lL+F+V18NIyNdVyJJkjRvGRYlLTybN7PP6GjXVUiSJM1rTkOVJEmSJA0wLEqSJEmSBhgWJUmSJEkDDIuSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNMCwKEmSJEkaYFiUJEmSJA0wLEqSJEmSBhgWJUmSJEkDDIuSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNMCwKEmSJEkakKrquoZZtWzVIbXq+FO7LmNGNqwd4+QtS7suQ5qzznzPiaxeUTz22Nd1XYo0p/l5Ik3N80Qzse2kY7ouYVJJLqqqw6fazpFFSZIkSdIAw6IkSZIkaYBhUZIkSZI0wLAoSZIkSRpgWJQkSZIkDTAsSpIkSZIGDDUsJlmZ5KwkVyTZmuRRSV6V5Jokm9ufo9tt1yT5n772f+47zqeSXJLk8iT/nGRJ37oXJvlqu+71w+yPJEmSJC0Ww/7SmDcBn6qqpyfZE9gbeDLwD1X1xkm2/2ZVHTZJ++9W1Q+TBDgLeAZwZpJfBY4FDq2q25LcfUj9kCRJkqRFZWhhMcldgCOBEwCq6ifAT5q8d8dU1Q/bh0uBPYFql/8EOKmqbmu3u35mVUuSJEmSYLjTUO8LfB94Z5IvJ3l7khXtuhckuTTJaUnu2rfPfdptz0vyuP6DJfk0cD1wC83oIsD9gcclubDd54gh9keSJEmSFo1hTkNdCjwMeGFVXZjkTcCJwFuAv6UZHfxb4GTgj4DrgHtX1Y1JHg58OMmDxkcVq+rJSfYCzgCOAs5tn+OuwCOBI4D3J7lvVVV/IUnWA+sBVu53ABvWjg2x28N34HLmfR+kYVq9othzieeJNBU/T6SpeZ5oJnq9XtclzMgww+LVwNVVdWG7fBZwYlV9b3yDJG8DPg7QTiUdn056UZJv0owcbhrfvqp+nOSjNNcpnts+x9ltOPxSkp8B+9OMaNK330ZgI8CyVYfUyVuGfanmcG1YO8Z874M0TEfcGlavKM8TaQp+nkhT8zzRTGxbN9J1CTMytGmoVfVd4Kokv9Q2PQH4SpJVfZs9DbgMIMkB43c5TXJf4BDgyiT7jO+TZClwNHBFu/+HaUYZSXJ/musZbxhWnyRJkiRpsRj2f5O8EDijvRPqlcCzgTcnOYxmGuo24LnttkcCr0kyBuwAnldVNyU5EPhokmXAEuA/gPGv1TgNOC3JZcBPgOMnTkGVJEmSJN1xQw2LVbUZOHxC83E72faDwAcnaf8ezfWIk+3zE+BZMyxTkiRJkjTBMO+GKkmSJEmapwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjTAsChJkiRJGpDF9rWEhx9+eG3atKnrMmak1+sxMjLSdRnS3DUywujoKCs3b+66EmlO8/NEmprniRaiJBdV1cSvOBzgyKIkSZIkaYBhUZIkSZI0wLAoSZIkSRpgWJQkSZIkDVh0N7hZtuqQWnX8qV2XMSMb1o5x8palXZchzVlnvudEVq8oHnvs67ouRZrT/DyRpuZ5MnzbTjqm6xIWHW9wI0mSJEn6hRkWJUmSJEkDDIuSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNMCwKEmSJEkaMLSwmGSvJF9KckmSy5O8um0/I8lXk1yW5LQke7Ttf55kc/tzWZIdSe7WrluZ5KwkVyTZmuRRbfthSS5o99mU5BHD6o8kSZIkLSbDHFm8DTiqqh4CHAY8JckjgTOAXwbWAsuB5wBU1Ruq6rCqOgx4BXBeVd3UHutNwKeq6peBhwBb2/bXA69u9/mbdlmSJEmSNENLh3Xgqipge7u4R/tTVfWJ8W2SfAlYPcnuzwTe225zF+BI4IT2uD8BfjL+NMBd2sf7Atfu1k5IkiRJ0iI1tLAIkGQJcBFwP+CfqurCvnV7AMcBL56wz97AU4AXtE33Bb4PvDPJQ9rjvbiqbgVeAnw6yRtpRkkfvZM61gPrAVbudwAb1o7ttj524cDlzPs+SMO0ekWx5xLPE2kqfp5IU/M8Gb5er9d1CdqJoYbFqtoBHJZkJfChJA+uqsva1W8Fzq+q/5qw228Cn++bgroUeBjwwqq6MMmbgBOBvwb+BPizqvpgkt8F3gE8cZI6NgIbAZatOqRO3jLUbg/dhrVjzPc+SMN0xK1h9YryPJGm4OeJNDXPk+Hbtm6k6xK0E7NyN9SqGgV6NCOGJHklcADw0kk2/33aKaitq4Gr+0Ylz6IJjwDHA2e3jz8AeIMbSZIkSdoNhnk31APaEUWSLKcZ8bsiyXOAJwPPrKqfTdhnX+DxwEfG26rqu8BVSX6pbXoC8JX28bXt9gBHAV8fUnckSZIkaVEZ5pj6KuD09rrFOwHvr6qPJxkDvg18MQnA2VX1mnafpwGfaa9H7PdC4IwkewJXAs9u2/8X8KYkS4Ef016XKEmSJEmamWHeDfVS4KGTtO/0OavqXcC7JmnfDBw+SfvngIfPpE5JkiRJ0qBZuWZRkiRJkjS/GBYlSZIkSQMMi5IkSZKkAYZFSZIkSdIAw6IkSZIkaUCqqusaZtXhhx9emzZt6rqMGen1eoyMjHRdhjR3jYwwOjrKys2bu65EmtP8PJGm5nmihSjJRVU18G0TEzmyKEmSJEkaYFiUJEmSJA0wLEqSJEmSBhgWJUmSJEkDFt0NbpatOqRWHX9q12XMyIa1Y5y8ZWnXZUhz1pnvOZHVK4rHHvu6rkuR5jQ/T6SpzbfzZNtJx3RdguYBb3AjSZIkSfqFGRYlSZIkSQMMi5IkSZKkAYZFSZIkSdIAw6IkSZIkaYBhUZIkSZI0wLAoSZIkSRow1LCY5LQk1ye5bEL7C5N8NcnlSV7ftj0iyeb255IkT+vb/s/abS9L8t4ke0043j8m2T7MvkiSJEnSYjLskcV3AU/pb0jyq8CxwKFV9SDgje2qy4DDq+qwdp9/SbI0yUHAi9p1DwaWAL/fd7zDgZVD7ockSZIkLSpDDYtVdT5w04TmPwFOqqrb2m2ub//9UVWNtdvsBVTfPkuB5UmWAnsD1wIkWQK8AXj50DohSZIkSYvQ0g6e8/7A45K8Fvgx8LKq+m+AJL8CnAYcDBzXhsdrkrwR+A7wP8Bnquoz7bFeAHy0qq5LstMnTLIeWA+wcr8D2LB2bKfbzgcHLmfe90EaptUrij2XeJ5IU/HzRJrafDtPer1e1yVoAekiLC4F7go8EjgCeH+S+1bjQuBBSR4AnJ7kk8Bymmmr9wFGgQ8keRbwH8AzgJGpnrCqNgIbAZatOqRO3tJFt3efDWvHmO99kIbpiFvD6hXleSJNwc8TaWrz7TzZtm6k6xK0gHTxzr8aOLuqCvhSkp8B+wPfH9+gqrYmuRV4ME1I/FZVfR8gydnAo4EfAPcDvtGOKu6d5BtVdb9Z7Y0kSZIkLUBdfHXGh4GjAJLcH9gTuCHJfdprEklyMPBLwDaa6aePTLJ3mlT4BGBrVZ1TVfeoqjVVtQb4kUFRkiRJknaPoY4sJnkvzTTR/ZNcDbyS5prE09qv0/gJcHxVVZLHAicm+SnwM+BPq+oGmiB5FnAxMAZ8mXZKqSRJkiRpOIYaFqvqmTtZ9axJtn038O6dHOeVNEFzV8+1zx0uUJIkSZI0qS6moUqSJEmS5jjDoiRJkiRpgGFRkiRJkjTAsChJkiRJGmBYlCRJkiQNSFV1XcOsOvzww2vTpk1dlzEjvV6PkZGRrsuQ5q6REUZHR1m5eXPXlUhzmp8n0tQ8T7QQJbmoqg6fajtHFiVJkiRJAwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjRgadcFzLYt19zMmhPP6bqMGdmwdowT5nkfpGE688obWb2iOMzzRNolP0+00Gw76ZiuS5AWFEcWJUmSJEkDDIuSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNMCwKEmSJEkaYFiUJEmSJA3oJCwm+bMklye5LMl7k+yV5AVJvpGkkuzft+1dk3woyaVJvpTkwW37vZL8Z5Kt7bFe3EVfJEmSJGkhmvWwmOQg4EXA4VX1YGAJ8PvA54EnAt+esMtfApur6lDgD4E3te1jwIaqegDwSOD5SR44C12QJEmSpAWvq2moS4HlSZYCewPXVtWXq2rbJNs+EPgsQFVdAaxJcmBVXVdVF7fttwBbgYNmpXpJkiRJWuCWzvYTVtU1Sd4IfAf4H+AzVfWZXexyCfDbwOeSPAI4GFgNfG98gyRrgIcCF052gCTrgfUAK/c7gA1rx2bekQ4duJx53wdpmFavKPZc4nkiTcXPEy00vV5vtx9z+/btQzmuNB/MelhMclfgWOA+wCjwgSTPqqp/28kuJwFvSrIZ2AJ8mWYK6vjx9gE+CLykqn442QGqaiOwEWDZqkPq5C2z3u3dasPaMeZ7H6RhOuLWsHpFeZ5IU/DzRAvNtnUju/2YvV6PkZHdf1xpPujiE+KJwLeq6vsASc4GHg1MGhbbAPjsdtsA32p/SLIHTVA8o6rOHn7pkiRJkrQ4dHHN4neARybZuw1/T6C53nBSSVYm2bNdfA5wflX9sN33HcDWqjpl6FVLkiRJ0iIy62Gxqi4EzgIupplWeidgY5IXJbma5nrES5O8vd3lAcDlSa4Afh0Y/4qMxwDHAUcl2dz+HD2bfZEkSZKkhaqTCxWq6pXAKyc0v7n9mbjtF4FDJmn/HJChFChJkiRJi1xXX50hSZIkSZrDDIuSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNKCTu6F2ae1B+7LppGO6LmNGer0e29aNdF2GNHdd8AZGR0fZNs/PdWnY/DyRJO2KI4uSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNGDR3eBmyzU3s+bEc7ouY0Y2rB3jhHneB2mYzrzyRlavKA7zPJF2yc+TucWbckmaaxxZlCRJkiQNMCxKkiRJkgYYFiVJkiRJAwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjSgs7CYZEmSLyf5eLt8nyQXJvl6kvcl2bNtPzjJZ5NcmqSXZHXfMe6d5DNJtib5SpI13fRGkiRJkhaWLkcWXwxs7Vt+HfAPVXUI8APgj9v2NwL/WlWHAq8B/r5vn38F3lBVDwAeAVw/9KolSZIkaRHoJCy2o4PHAG9vlwMcBZzVbnI68Fvt4wcCn20f/ydwbLvPA4GlVXUuQFVtr6ofzUoHJEmSJGmBW9rR854KvBy4c7u8HzBaVWPt8tXAQe3jS4DfAd4EPA24c5L9gPsDo0nOBu4D/DtwYlXtmPhkSdYD6wFW7ncAG9aOTdxkXjlwOfO+D9IwrV5R7LnE80Saip8nc0uv1+u6BE1i+/bt/m60aM16WEzyG8D1VXVRkpHx5kk2rfbflwFvSXICcD5wDTBGU/vjgIcC3wHeB5wAvGPgQFUbgY0Ay1YdUidv6Soj7x4b1o4x3/sgDdMRt4bVK8rzRJqCnydzy7Z1I12XoEn0ej1GRka6LkPqRBefEI8BnprkaGAv4C40I40rkyxtRxdXA9cCVNW1wG8DJNkH+J2qujnJ1cCXq+rKdt2HgUcySViUJEmSJN0xs37NYlW9oqpWV9Ua4PeB/6iqdTTXIz693ex44CMASfZPMl7nK4DT2sf/Ddw1yQHt8lHAV2ahC5IkSZK04M2l71n8C+ClSb5Bcw3j+AjhCPDVJF8DDgReC9Bem/gy4LNJttBMZX3bbBctSZIkSQtRpxcqVFUP6LWPr6T5+ouJ25zFz++SOnHducChw6tQkiRJkhanuTSyKEmSJEmaIwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjSg07uhdmHtQfuy6aRjui5jRnq9HtvWjXRdhjR3XfAGRkdH2TbPz3Vp2Pw8kSTtiiOLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpgGFRkiRJkjTAsChJkiRJGrDo7oa65ZqbWXPiOV2XMSMb1o5xwjzvgzRMZ155I6tXFId5nki7NNc+T7yDsSTNLY4sSpIkSZIGGBYlSZIkSQMMi5IkSZKkAYZFSZIkSdIAw6IkSZIkaYBhUZIkSZI0wLAoSZIkSRrQSVhMcq8k/5lka5LLk7y4bb9bknOTfL39964T9jsiyY4kT+9re317jK1J3pwks90fSZIkSVpouhpZHAM2VNUDgEcCz0/yQOBE4LNVdQjw2XYZgCRLgNcBn+5rezTwGOBQ4MHAEcDjZ6sTkiRJkrRQdRIWq+q6qrq4fXwLsBU4CDgWOL3d7HTgt/p2eyHwQeD6/kMBewF7AsuAPYDvDbV4SZIkSVoElnZdQJI1wEOBC4EDq+o6aAJlkru32xwEPA04imb0kHabLyb5T+A6IMBbqmrrJM+xHlgPsHK/A9iwdmyYXRq6A5cz7/sgDdPqFcWeSzxPpKnMtc+TXq/XdQnSgO3bt/ve1KLVaVhMsg/NaOFLquqHu7jc8FTgL6pqR/82Se4HPABY3Tadm+TIqjq/f+eq2ghsBFi26pA6eUvnGXlGNqwdY773QRqmI24Nq1eU54k0hbn2ebJt3UjXJUgDer0eIyMjXZchdaKzT4gke9AExTOq6uy2+XtJVrWjiqv4+ZTTw4Ez26C4P3B0kjHgEOCCqtreHvOTNNdA3i4sSpIkSZLumK7uhhrgHcDWqjqlb9VHgePbx8cDHwGoqvtU1ZqqWgOcBfxpVX0Y+A7w+CRL2/D5eJrrHyVJkiRJM9DVyOJjgOOALUk2t21/CZwEvD/JH9MEwWdMcZyzaK5j3EJzs5tPVdXHhlOyJEmSJC0enYTFqvoczQ1pJvOEKfY9oe/xDuAfOWZxAAAgAElEQVS5u68ySZIkSRJ09z2LkiRJkqQ5zLAoSZIkSRpgWJQkSZIkDTAsSpIkSZIGGBYlSZIkSQO6+uqMzqw9aF82nXRM12XMSK/XY9u6ka7LkOauC97A6Ogo2+b5uS4Nm58nkqRdcWRRkiRJkjTAsChJkiRJGmBYlCRJkiQNMCxKkiRJkgYsuhvcbLnmZtaceE7XZczIhrVjnDDP+yAN05lX3sjqFcVhnieap7w5kyRpLnBkUZIkSZI0wLAoSZIkSRpgWJQkSZIkDTAsSpIkSZIGGBYlSZIkSQMMi5IkSZKkAYZFSZIkSdKATsJiktOSXJ/ksr62uyU5N8nX23/v2ravS3Jp+/OFJA+ZcKwlSb6c5OOz3Q9JkiRJWqi6Gll8F/CUCW0nAp+tqkOAz7bLAN8CHl9VhwJ/C2ycsN+Lga3DK1WSJEmSFp9OwmJVnQ/cNKH5WOD09vHpwG+1236hqn7Qtl8ArB7fIclq4Bjg7UMtWJIkSZIWmaVdF9DnwKq6DqCqrkty90m2+WPgk33LpwIvB+68qwMnWQ+sB1i53wFsWDu2eyruyIHLmfd9kIZp9YpizyWeJ5q/er3erDzP9u3bZ+25pPnK80SL2VwKi7uU5FdpwuJj2+XfAK6vqouSjOxq36raSDt9ddmqQ+rkLfOm25PasHaM+d4HaZiOuDWsXlGeJ5q3tq0bmZXn6fV6jIzMznNJ85XniRazuXQ31O8lWQXQ/nv9+Iokh9JMNT22qm5smx8DPDXJNuBM4Kgk/za7JUuSJEnSwjSXwuJHgePbx8cDHwFIcm/gbOC4qvra+MZV9YqqWl1Va4DfB/6jqp41uyVLkiRJ0sLUyRytJO8FRoD9k1wNvBI4CXh/kj8GvgM8o938b4D9gLcmARirqsNnvWhJkiRJWkQ6CYtV9cydrHrCJNs+B3jOFMfrAb0ZFyZJkiRJAubWNFRJkiRJ0hxhWJQkSZIkDTAsSpIkSZIGGBYlSZIkSQMMi5IkSZKkAZ3cDbVLaw/al00nHdN1GTPS6/XYtm6k6zKkueuCNzA6Osq2eX6uS5IkdcmRRUmSJEnSAMOiJEmSJGnAHQqLSe6U5C7DKkaSJEmSNDdMGRaTvCfJXZKsAL4CfDXJnw+/NEmSJElSV6YzsvjAqvoh8FvAJ4B7A8cNtSpJkiRJUqemczfUPZLsQRMW31JVP01SQ65raLZcczNrTjyn6zJmZMPaMU6Y532QhunMK29k9YriMM+TBc873kqSNDzTGVn8F2AbsAI4P8nBwA+HWZQkSZIkqVtTjixW1ZuBN/c1fTvJrw6vJEmSJElS13YaFpO8dIp9T9nNtUiSJEmS5ohdjSzeedaqkCRJkiTNKTsNi1X16tksRJIkSZI0d0x5zWKSvYA/Bh4E7DXeXlV/NMS6JEmSJEkdms7dUN8N3AN4MnAesBq4ZZhFSZIkSZK6NZ2weL+q+mvg1qo6HTgGWDusgpJsS7IlyeYkm9q2ZyS5PMnPkhzet+2vJbmo3f6iJEcNqy5JkiRJWkymnIYK/LT9dzTJg4HvAmuGVlHjV6vqhr7ly4DfpvnOx343AL9ZVde2tX0aOGjItUmSJEnSgjedsLgxyV2BvwY+CuwD/M1Qq5qgqrYCJJnY/uW+xcuBvZIsq6rbZrE8SZIkSVpwpgyLVfX29uF5wH2HW07zlMBnkhTwL1W1cZr7/Q7w5cmCYpL1wHqAlfsdwIa1Y7ut2C4cuJx53wdpmFavKPZc4nmyGPR6va5LmNe2b9/uayhNwfNEi9l07oa6jCaIrenfvqpeM6SaHtNOK707cG6SK6rq/ClqfBDwOuBJk61vA+dGgGWrDqmTt0xnQHXu2rB2jPneB2mYjrg1rF5RnieLwLZ1I12XMK/1ej1GRka6LkOa0zxPtJhN5wY3HwGOBcaAW/t+hqKqrm3/vR74EPCIXW2fZHW73R9W1TeHVZckSZIkLSbT+W/31VX1lKFXAiRZAdypqm5pHz8J2OkIZpKVwDnAK6rq87NRoyRJkiQtBtMZWfxCkqF9VcYEBwKfS3IJ8CXgnKr6VJKnJbkaeBRwTpJPt9u/ALgf8NftV21sbqevSpIkSZJmYDoji48FTkjyLeA2IEBV1aG7u5iquhJ4yCTtH6KZajqx/e+Av9vddUiSJEnSYjedsPjrQ69CkiRJkjSnTCcs3jLNNkmSJEnSAjGdaxYvBr4PfA34evv4W0kuTvLwYRYnSZIkSerGdMLip4Cjq2r/qtqPZlrq+4E/Bd46zOIkSZIkSd2YTlg8vKrG7z5KVX0GOLKqLgCWDa0ySZIkSVJnpnPN4k1J/gI4s13+PeAHSZYAPxtaZUOy9qB92XTSMV2XMSO9Xo9t60a6LkOauy54A6Ojo2yb5+e6JElSl6YzsvgHwGrgw8BHgHu3bUuA3x1eaZIkSZKkrkw5slhVNwAv3Mnqb+zeciRJkiRJc8FOw2KSU6vqJUk+BtTE9VX11KFWJkmSJEnqzK5GFt/d/vvG2ShEkiRJkjR37DQsVtVF7b/njbcluStwr6q6dBZqG4ot19zMmhPP6bqMGdmwdowT5nkfpOnwBjWSJEndmfIGN0l6Se6S5G7AJcA7k5wy/NIkSZIkSV2Zzt1Q962qHwK/Dbyzqh4OPHG4ZUmSJEmSujSdsLg0ySqar8n4+JDrkSRJkiTNAdMJi68BPg18o6r+O8l9ga8PtyxJkiRJUpem8z2LHwA+0Ld8JfA7wyxKkiRJktSt6YwsSpIkSZIWGcOiJEmSJGnATsNikhe3/z5mtopJcq8k/5lka5LL+2p4VZJrkmxuf47u2+fQJF9st9+SZK/ZqleSJEmSFqpdXbP4bOBNwD8CD5udchgDNlTVxUnuDFyU5Nx23T9U1Rv7N06yFPg34LiquiTJfsBPZ6lWSZIkSVqwdhUWtybZBhyQ5NK+9gBVVYfu7mKq6jrguvbxLUm2AgftYpcnAZdW1SXtPjfu7pokSZIkaTHaaVisqmcmuQfN12Y8dfZKaiRZAzwUuBB4DPCCJH8IbKIZffwBcH+gknwaOAA4s6peP8mx1gPrAVbudwAb1o7NSh+G5cDlzPs+SNPR6/V+of0OGx1lx44dv/D+0mKxfft2zxNpCp4nWsxSVVNvlOxJE8wAvlpVQ53qmWQf4DzgtVV1dpIDgRuAAv4WWFVVf5TkZcDzgSOAHwGfBf6qqj67s2MvW3VIrTr+1GGWP3Qb1o5x8pYpv/VEmve2nXTML7bjyAijo6Os3Lx59xYkLTC9Xo+RkZGuy5DmNM8TLURJLqqqw6fabsq7oSZ5PPB14J+AtwJfS3LkzEvc6fPtAXwQOKOqzgaoqu9V1Y6q+hnwNuAR7eZXA+dV1Q1V9SPgE8ze9ZWSJEmStGBN56szTgGeVFWPr6ojgScD/zCMYpIEeAewtapO6Wtf1bfZ04DL2sefBg5Nsnd7s5vHA18ZRm2SJEmStJhMZy7jHlX11fGFqvpaO/o3DI8BjgO2JBmfP/aXwDOTHEYzDXUb8Ny2lh8kOQX473bdJ6rqnCHVJkmSJEmLxnTC4qYk7wDe3S6vAy4aRjFV9Tmau61O9Ild7PNvNF+fIUmSJEnaTaYTFv+E5iYyL6IJcufTXLsoSZIkSVqgpgyLVXUbzXWLp0y1rSRJkiRpYZjODW4kSZIkSYuMYVGSJEmSNMCwKEmSJEkaMOU1i0nuD/w5cHD/9lV11BDrGpq1B+3LppOO6bqMGen1emxbN9J1GZIkSZIWsOncDfUDwD8DbwN2DLccSZIkSdJcMJ2wOFZV///QK5EkSZIkzRk7DYtJ7tY+/FiSPwU+BNw2vr6qbhpybZIkSZKkjuxqZPEioIC0y3/et66A+w6rKEmSJElSt3YaFqvqPgBJ9qqqH/evS7LXsAuTJEmSJHVnOtcsfgF42DTa5oUt19zMmhPP6bqMGdmwdowT5nkftLBsm+d3GJYkSdKgXV2zeA/gIGB5kofy8+modwH2noXaJEmSJEkd2dXI4pOBE4DVwCl97bcAfznEmiRJkiRJHdvVNYunA6cn+Z2q+uAs1iRJkiRJ6th0rlk8OMlLJ7TdDFxUVZuHUJMkSZIkqWN3msY2hwPPo7l+8SBgPTACvC3Jy4dXmiRJkiSpK9MZWdwPeFhVbQdI8krgLOBImu9ifP3wypMkSZIkdWE6I4v3Bn7St/xT4OCq+h/gtt1ZTJLTklyf5LK+tjckuSLJpUk+lGRl275HktOTbEmyNckrdmctkiRJkrSYTScsvge4IMkr21HFzwPvTbIC+MpuruddwFMmtJ0LPLiqDgW+BoyHwmcAy6pqLfBw4LlJ1uzmeiRJkiRpUZpyGmpV/W2STwKPofmuxedV1aZ29brdWUxVnT8x8FXVZ/oWLwCePr4KWJFkKbCcZvTzh7uzHkmSJElarKZzzSLAl4Frx7dPcu+q+s7Qqtq5PwLe1z4+CzgWuA7YG/izqrqpg5okSZIkacGZMiwmeSHwSuB7wA6a0cUCDh1uaQN1/G9gDDijbXpEW889gbsC/5Xk36vqykn2XU9zF1dW7ncAG9aOzU7RQ3LgcuZ9H7Sw9Hq9rku4ncNGR9mxY8ecq0uaa7Zv3+55Ik3B80SL2XRGFl8M/FJV3TjsYnYmyfHAbwBPqKpqm/8A+FRV/RS4Psnnab7mYyAsVtVGYCPAslWH1MlbpjugOjdtWDvGfO+DFpZt60a6LuH2Vq5kdHSUkZGRriuR5rRer+d5Ik3B80SL2XRucHMVcPOwC9mZJE8B/gJ4alX9qG/Vd4Cj0lgBPBK4oosaJUmSJGmhmc7w1JVAL8k59H1VRlWdsruLSfJeYATYP8nVNNNfXwEsA85NAnBBVT0P+CfgncBlNFNj31lVl+7umiRJkiRpMZpOWPxO+7Nn+zM0VfXMSZrfsZNtt9N8fYYkSZIkaTebzldnvBogyYqqunX4JUmSJEmSujblNYtJHpXkK8DWdvkhSd469MokSZIkSZ2Zzg1uTgWeDNwIUFWXAEcOsyhJkiRJUremExapqqsmNO0YQi2SJEmSpDliOje4uSrJo4FKsifwItopqZIkSZKkhWk6I4vPA54PHARcDRwG/Okwi5IkSZIkdWs6d0O9AVjX35bkJTTXMs47aw/al00nHdN1GTPS6/XYtm6k6zIkSZIkLWDTumZxEi/drVVIkiRJkuaUXzQsZrdWIUmSJEmaU37RsFi7tQpJkiRJ0pyy02sWk9zC5KEwwPKhVSRJkiRJ6txOw2JV3Xk2C5ktW665mTUnntN1GTOyYe0YJ8zzPmhu2TbPb/okSZKk3e8XnYYqSZIkSVrADIuSJEmSpAGGRUmSJEnSAMOiJEmSJGmAYVGSJEmSNMCwKEmSJEkaYFiUJEmSJA2YN2ExycokZyW5IsnWJI/qW/eyJJVk/y5rlCRJkqSFYmnXBdwBbwI+VVVPT7InsDdAknsBvwZ8p8viJEmSJGkhmRcji0nuAhwJvAOgqn5SVaPt6n8AXg5UR+VJkiRJ0oIzX0YW7wt8H3hnkocAFwEvBp4AXFNVlyTZ6c5J1gPrAVbudwAb1o4Nv+IhOnA5874Pmlt6vV7XJexWh42OsmPHjgXXL2l32759u+eJNAXPEy1m8yUsLgUeBrywqi5M8ibgVTSjjU+aaueq2ghsBFi26pA6ect86fbkNqwdY773QXPLtnUjXZewe61cyejoKCMjI11XIs1pvV7P80SagueJFrN5MQ0VuBq4uqoubJfPogmP9wEuSbINWA1cnOQe3ZQoSZIkSQvHvAiLVfVd4Kokv9Q2PQG4uKruXlVrqmoNTaB8WLutJEmSJGkG5tNcxhcCZ7R3Qr0SeHbH9UiSJEnSgjVvwmJVbQYO38X6NbNXjSRJkiQtbPNiGqokSZIkaXYZFiVJkiRJAwyLkiRJkqQBhkVJkiRJ0gDDoiRJkiRpwLy5G+rusvagfdl00jFdlzEjvV6PbetGui5DkiRJ0gLmyKIkSZIkaYBhUZIkSZI0wLAoSZIkSRpgWJQkSZIkDTAsSpIkSZIGLLq7oW655mbWnHhO12XMyIa1Y5wwz/ugO2bbPL+DryRJkuYfRxYlSZIkSQMMi5IkSZKkAYZFSZIkSdIAw6IkSZIkaYBhUZIkSZI0wLAoSZIkSRowL8Jikr2SfCnJJUkuT/Lqtv2MJF9NclmS05Ls0XWtkiRJkrQQzIuwCNwGHFVVDwEOA56S5JHAGcAvA2uB5cBzuitRkiRJkhaOpV0XMB1VVcD2dnGP9qeq6hPj2yT5ErC6g/IkSZIkacGZLyOLJFmSZDNwPXBuVV3Yt24P4DjgU13VJ0mSJEkLybwYWQSoqh3AYUlWAh9K8uCquqxd/Vbg/Kr6r8n2TbIeWA+wcr8D2LB2bFZqHpYDlzPv+6A7ptfrdV3CvHLY6Cg7duzwdZOmsH37ds8TaQqeJ1rM5k1YHFdVo0l6wFOAy5K8EjgAeO4u9tkIbARYtuqQOnnLvOv27WxYO8Z874PumG3rRrouYX5ZuZLR0VFGRka6rkSa03q9nueJNAXPEy1m82IaapID2hFFkiwHnghckeQ5wJOBZ1bVz7qsUZIkSZIWkvkyPLUKOD3JEpqA+/6q+niSMeDbwBeTAJxdVa/psE5JkiRJWhDmRVisqkuBh07SPi/qlyRJkqT5Zl5MQ5UkSZIkzS7DoiRJkiRpgGFRkiRJkjTAsChJkiRJGmBYlCRJkiQNMCxK0v9t7/6DNbvr+oC/P5OVJLR0V0PMpBtlSQk4wLWLpY4VsZeAbSQKtjACs1BDaXeIViizTo2tI850rDuUKMZqmbXGBYREiRlEooUWfUwVFEJANyFCmbBogtOIzE1dh0B3+frHfVI297u7997s3Xue8zyv18ydufc858f7PLsfnn1zzj0BAKCzcP/piaXdO3PHwauHjnFWJpNJju5bHjoGAAAwx1xZBAAAoKMsAgAA0FEWAQAA6CiLAAAAdBbuATdH7n8we667begYZ+XA0vFcM/JzWCRHR/5AJQAAFpMriwAAAHSURQAAADrKIgAAAB1lEQAAgI6yCAAAQEdZBAAAoKMsAgAA0BlFWayqG6vqgaq6a83yH6yqT1TV3VX1hqHyAQAAzJtRlMUkh5NcdfKCqnpOkhcm+cbW2tOSvHGAXAAAAHNpFGWxtXZ7ks+vWXxtkoOttS9O13lg24MBAADMqR1DBzgLT07y7Kr6iSQPJfmh1tqHT7ViVe1Psj9Jdl10cQ4sHd++lOfAJRdm9OewSCaTydARFs7elZWcOHHCew/rOHbsmDmBdZgTFtmYy+KOJF+d5FuS/MMkv1pVl7fW2toVW2uHkhxKkvMvvaJdf2TMp71aFMd+Dovk6L7loSMsnl27srKykuXl5aGTwEybTCbmBNZhTlhko7gN9TTuS3JrW/WhJF9O8viBMwEAAMyFMZfFdyW5Mkmq6slJHpPkc4MmAgAAmBOjuJexqm5Kspzk8VV1X5LXJ7kxyY3T/5zGl5J836luQQUAAGDzRlEWW2svO81LL9/WIAAAAAtizLehAgAAcI4oiwAAAHSURQAAADrKIgAAAB1lEQAAgM4onoa6lZZ278wdB68eOsZZmUwmObpveegYAADAHHNlEQAAgI6yCAAAQEdZBAAAoKMsAgAA0Fm4B9wcuf/B7LnutqFjnJUDS8dzzcjPYV4cHfnDkgAA4HRcWQQAAKCjLAIAANBRFgEAAOgoiwAAAHSURQAAADrKIgAAAB1lEQAAgM7oy2JVva6q7q6qu6rqpqq6YOhMAAAAYzfqslhVu5O8JskzW2tPT3JekpcOmwoAAGD8Rl0Wp3YkubCqdiR5bJLPDpwHAABg9HYMHeBstNbur6o3JvnTJF9I8r7W2vvWrldV+5PsT5JdF12cA0vHtzfoFrvkwoz+HObFZDIZOgKnsHdlJSdOnPDnA+s4duyYOYF1mBMW2ajLYlV9dZIXJnlikpUk76yql7fWfvnk9Vprh5IcSpLzL72iXX9k1KedA0vHM/ZzmBdH9y0PHYFT2bUrKysrWV5eHjoJzLTJZGJOYB3mhEU29ttQn5fk0621v2it/b8ktyb51oEzAQAAjN7Yy+KfJvmWqnpsVVWS5ya5Z+BMAAAAozfqstha+8MktyS5M8mRrJ7PoUFDAQAAzIHR/+Jba+31SV4/dA4AAIB5MuoriwAAAJwbyiIAAAAdZREAAICOsggAAEBHWQQAAKAz+qehbtbS7p254+DVQ8c4K5PJJEf3LQ8dAwAAmGOuLAIAANBRFgEAAOgoiwAAAHSURQAAADrKIgAAAJ2FexrqkfsfzJ7rbhs6xlk5sHQ814z8HMbu6MifqAsAAOtxZREAAICOsggAAEBHWQQAAKCjLAIAANBRFgEAAOgoiwAAAHSURQAAADpzURar6ryq+mhVvWfoLAAAAPNgLspiktcmuWfoEAAAAPNi9GWxqi5LcnWS/zZ0FgAAgHmxY+gAW+BNSf5dksedboWq2p9kf5LsuujiHFg6vk3Rzo1LLszoz2HsJpPJ0BE4g70rKzlx4oQ/J1jHsWPHzAmsw5ywyEZdFqvqu5I80Fr7SFUtn2691tqhJIeS5PxLr2jXHxn1aefA0vGM/RzG7ui+5aEjcCa7dmVlZSXLy8tDJ4GZNplMzAmsw5ywyMZ+G+qzkrygqo4muTnJlVX1y8NGAgAAGL9Rl8XW2o+01i5rre1J8tIkv91ae/nAsQAAAEZv1GURAACAc2NufvGttTZJMhk4BgAAwFxwZREAAICOsggAAEBHWQQAAKCjLAIAANBRFgEAAOjMzdNQN2pp987ccfDqoWOclclkkqP7loeOAQAAzDFXFgEAAOgoiwAAAHSURQAAADrKIgAAAJ2Fe8DNkfsfzJ7rbhs6xlk5sHQ814z8HNhaR0f+0CYAAGaPK4sAAAB0lEUAAAA6yiIAAAAdZREAAICOsggAAEBHWQQAAKCjLAIAANAZdVmsqq+rqt+pqnuq6u6qeu3QmQAAAObBjqEDnKXjSQ601u6sqscl+UhV/Y/W2seHDgYAADBmo76y2Fr789bandPv/yrJPUl2D5sKAABg/MZ+ZfH/q6o9SZ6R5A9P8dr+JPuTZNdFF+fA0vFtzbbVLrkwoz8HttZkMhk6wkzZu7KSEydOeF9gHceOHTMnsA5zwiKbi7JYVX87ya8l+bettf+79vXW2qEkh5Lk/EuvaNcfGfdpH1g6nrGfA1vr6L7loSPMll27srKykuXl5aGTwEybTCbmBNZhTlhko74NNUmq6quyWhTf3lq7deg8AAAA82DUZbGqKskvJrmntfZTQ+cBAACYF6Mui0meleQVSa6sqo9Nv54/dCgAAICxG/UvvrXWfi9JDZ0DAABg3oz9yiIAAADngLIIAABAR1kEAACgoywCAADQURYBAADojPppqI/G0u6duePg1UPHOCuTySRH9y0PHQMAAJhjriwCAADQURYBAADoKIsAAAB0lEUAAAA6yiIAAACdhXsa6pH7H8ye624bOkbn6Mif0AoAAMwXVxYBAADoKIsAAAB0lEUAAAA6yiIAAAAdZREAAICOsggAAEBHWQQAAKAz+rJYVTdW1QNVddfQWQAAAObF6MtiksNJrho6BAAAwDwZfVlsrd2e5PND5wAAAJgnO4YOsB2qan+S/Umy66KLc2Dp+MCJepPJZMPrHjt2bFPrw6LZu7KSEydOmBNYh88TWJ85YZEtRFlsrR1KcihJzr/0inb9kdk77aP7lje87mQyyfLyxteHhbNrV1ZWVswJrMPnCazPnLDIRn8bKgAAAFtPWQQAAKAz+rJYVTcl+WCSp1TVfVX1qqEzAQAAjN3s/fLeJrXWXjZ0BgAAgHkz+iuLAAAAbD1lEQAAgI6yCAAAQEdZBAAAoKMsAgAA0Bn901A3a2n3ztxx8OqhYwAAAMw0VxYBAADoKIsAAAB0lEUAAAA6yiIAAAAdZREAAICOsggAAEBHWQQAAKCjLAIAANBRFgEAAOgoiwAAAHSURQAAADrKIgAAAB1lEQAAgI6yCAAAQEdZBAAAoKMsAvNn794ce9KThk4BADBqO4YOALDl3vSmfGoyyWVD5wAAGDFXFgEAAOgoiwAAAHSURQAAADrKIgAAAB1lEQAAgI6yCAAAQEdZBAAAoKMsAgAA0FEWAQAA6CiLAAAAdJRFAAAAOsoiAAAAHWURAACAjrIIAABAR1kEAACgoywCAADQURYBAADoKIsAAAB0lEUAAAA6yiIAAAAdZREAAICOsggAAEBHWQQAAKBTrbWhM2yrqvqLJJ85adHOJA9ucjeb3WYz629k3ccn+dwmjj8PHs2f07myHVm2+hhnu79Hu/1W/93fzLrmZFjblWUrj2NOFoM5GXZf5mQczMmw+9qOOXlCa+3idddqrS30V5JD53qbzay/kXWT3DH0+zaGP6cxZ9nqY5zt/h7t9lv9d38z65qTxciylccxJ4vxZU6G3Zc5GceXORl2X9sxJxv9chtq8hvbsM1m1n80eRbBLL0v25Flq49xtvt7tNufq7/7s/T3YZbM0vuyXVm28jjmZDHM0vtiTs7Ndubk7M3S+2JOtu+4nYW7DXUeVNUdrbVnDp0DZpk5gfWZE1ifOWGRubI4ToeGDgAjYE5gfeYE1mdOWFiuLAIAANBxZREAAICOsggAAEBHWQQAAKCjLAIAANBRFudMVV1eVb9YVbcMnQVmRVX9rap6S1X9QlXtGzoPzCqfIbC+qvqe6efJr1fVPxk6D5xLyuIMqaobq+qBqrprzfKrquoTVfWpqrruTPtord3bWnvVuU0Kw9vkvPzzJLe01v51khdse1gY0GZmxWcIi2qTc/Ku6efJNUleMkBc2DbK4mw5nOSqkxdU1XlJfi7JdyZ5apKXVdVTq2qpqt6z5utrtz8yDOZwNjgvSS5L8mfT1U5sY0aYBYez8VmBRXU4m5+TH52+Dv777lAAAAZBSURBVHNrx9AB+IrW2u1VtWfN4m9O8qnW2r1JUlU3J3lha+0nk3zX9iaE2bGZeUlyX1YL48fi/yRjwWxyVj6+velgNmxmTqrqniQHk/xWa+3ObQ0K28w/mmbf7nzlikiy+o/e3adbuaouqqo3J3lGVf3IuQ4HM+Z083JrkhdV1X9N8htDBIMZc8pZ8RkCj3C6z5QfTPK8JC+uqlcPEQy2iyuLs69OsaydbuXW2l8m8T9cLKpTzktr7a+TvHK7w8AMO92s+AyBrzjdnNyQ5IbtDgNDcGVx9t2X5OtO+vmyJJ8dKAvMOvMCG2NWYH3mhIWnLM6+Dye5oqqeWFWPSfLSJO8eOBPMKvMCG2NWYH3mhIWnLM6QqropyQeTPKWq7quqV7XWjif5N0nem+SeJL/aWrt7yJwwC8wLbIxZgfWZEzi1au20v/4GAADAgnJlEQAAgI6yCAAAQEdZBAAAoKMsAgAA0FEWAQAA6CiLAAAAdJRFAGZCVV1SVe+oqnur6iNV9cGq+mfrbPN3q+qWLTr+s6vq7qr6WFVduOa111TVPVX19k3uc1dVff9W5NsqVbVcVe9ZZ51H5N7K9xmA8VAWARhcVVWSdyW5vbV2eWvtHyR5aZLLzrRda+2zrbUXb1GMfUne2Frb21r7wprXvj/J81tr+za5z13TbcfmEbm3+H0GYCSURQBmwZVJvtRae/PDC1prn2mt/WySVNWeqvpfVXXn9OtbT1p+1/T7a6rq1qr671X1v6vqDac6UFU9t6o+WlVHqurGqjq/qv5Vku9N8mNrrx5W1ZuTXJ7k3VX1uqr65qr6wHQfH6iqp0zXe1pVfWh6ZfKPq+qKJAeT/L3psv98iiz/YrruH1XV26bLnlBV758uf39Vff10+eGqumF6zHur6sXT5b9SVc8/aZ+Hq+pFVXVBVf3S9Dw/WlXPOcXxf7yqfuikn++qqj1rc695n0+5342+/wCMx46hAwBAkqclufMMrz+Q5Dtaaw9NS9hNSZ55ivX2JnlGki8m+URV/Wxr7c8efrGqLkhyOMlzW2ufrKq3Jrm2tfamqvq2JO9prT3idsvW2qur6qokz2mtfa6q/k6Sb2+tHa+q5yX5T0lelOTVSX6mtfb2qnpMkvOSXJfk6a21vWuDVtXTkvyHJM+a7vdrpi/9lyRvba29par+ZZIbknzP9LVLk3xbkm9I8u4ktyS5OclLkvzm9LjPTXJtkh+Y5l+qqm9I8r6qevIZ3uOTPSL3tEA+7Ez7PeP7D8C4uLIIwMypqp+bXm378HTRVyX5hao6kuSdSZ56mk3f31p7sLX2UJKPJ3nCmtefkuTTrbVPTn9+S5Jv32S8nUneOb3S9tNZLbpJ8sEk/76qfjjJE05xK+taVya5pbX2uSRprX1+uvwfJXnH9Pu3ZbUcPuxdrbUvt9Y+nuSS6bLfSnJlVZ2f5DuzeivvF6bbvW267z9J8pkkGy2LZ3Km/a73/gMwIsoiALPg7iTf9PAPrbUfyOoVsouni16X5P8k+ftZvaL4mNPs54snfX8i/R00tQVZ/2OS32mtPT3Jdye5YJr5HUlekOQLSd5bVVeus59K0jZwvJPXOfn8anrch5JMkvzTrF5hvPnk19dxPI/8t8AFG9jmTPtd7/0HYESURQBmwW8nuaCqrj1p2WNP+n5nkj9vrX05ySuyeovno/EnSfZU1ZOmP78iye9uch87k9w//f6ahxdW1eVJ7m2t3ZDVW0S/MclfJXncafbz/iTfW1UXTbd/+DbUD2T14T7J6kN3fm8DmW5O8sokz07y3umy26fbZ3qb6Ncn+cSa7Y5mWtKr6puSPHG6/Ey5N7JfAOaAsgjA4FprLau/l/ePq+rTVfWhrN4i+sPTVX4+yfdV1R9k9ZbHv36Ux3koq6XqndNbWr+c5M1n3qrzhiQ/WVW/n0eW1pckuauqPpbV3yl8a2vtL5P8/vTBMY94wE1r7e4kP5Hkd6vqj5L81PSl1yR5ZVX9cVbL7Gs3kOl9Wb2d9n+21r40XfbzSc6bnuevJLmmtfbFNdv9WpKvmWa+Nsknp9lOm3uD+wVgDtTq5zMAAAB8hSuLAAAAdJRFAAAAOsoiAAAAHWURAACAjrIIAABAR1kEAACgoywCAADQ+RsQxkL0H+tEtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import timeit\n",
    "\n",
    "n = np.arange(17)  # lengths = 2**n to evaluate\n",
    "reps = 50  # number of repetitions for timeit\n",
    "\n",
    "gain = np.zeros(len(n))\n",
    "for N in n:\n",
    "    length = 2**N\n",
    "    # setup environment for timeit\n",
    "    tsetup = 'import numpy as np; from numpy.fft import rfft, irfft; \\\n",
    "            x=np.random.randn(%d); h=np.random.randn(%d)' % (length, length)\n",
    "    # direct convolution\n",
    "    tc = timeit.timeit('np.convolve(x, x, mode=\"full\")', setup=tsetup, number=reps)\n",
    "    # fast convolution\n",
    "    tf = timeit.timeit('irfft(rfft(x, %d) * rfft(h, %d))' % (2*length, 2*length), setup=tsetup, number=reps)\n",
    "    # speedup by using the fast convolution\n",
    "    gain[N] = tc/tf\n",
    "\n",
    "# show the results\n",
    "plt.figure(figsize = (15, 10))\n",
    "plt.barh(n, gain, log=True)\n",
    "plt.plot([1, 1], [-1, n[-1]+1], 'r-')\n",
    "plt.yticks(n, 2**n)\n",
    "plt.xlabel('Gain of fast convolution')\n",
    "plt.ylabel('Length of signals')\n",
    "plt.title('Comparison between direct/fast convolution')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* When is the fast convolution more efficient/faster than a direct convolution? \n",
    "* Why is it slower below a given signal length?\n",
    "* Is the trend of the gain as expected by the numerical complexity of the FFT?\n",
    "\n",
    "Solution: The gain in execution time of a fast convolution over a direct implementation of the the convolution for different signal lengths depends heavily on the particular implementation and hardware used. The fast convolution in this example is faster for two signals having a length equal or larger than 1024 samples. Discarding the outliers and short lengths, the overall trend in the gain is approximately logarithmic as predicted above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "**Copyright**\n",
    "\n",
    "This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2018*."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
